{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "\n",
    "from seqm.seqm_functions.constants import Constants\n",
    "from seqm.basics import Energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time : 7.856040 sec\n"
     ]
    }
   ],
   "source": [
    "#check code to produce energy terms for each molecule\n",
    "# with a 'learned' given parameters\n",
    "\n",
    "\n",
    "torch.set_default_dtype(torch.float64)\n",
    "if torch.cuda.is_available():\n",
    "    device = torch.device('cuda')\n",
    "else:\n",
    "    device = torch.device('cpu')\n",
    "\n",
    "N = 100\n",
    "at = 0\n",
    "#dir = torch.tensor([1.0,0.0,1.0],dtype=dtype,device=device)\n",
    "#dir=torch.randn(3,dtype=dtype,device=device)\n",
    "#dir /= torch.norm(dir)\n",
    "#print(dir.tolist())\n",
    "dir = torch.tensor([-0.3651281073742882, -0.5962853405673783, 0.7149302468281195],device=device)\n",
    "dxmin = -0.5\n",
    "dxmax = 0.5\n",
    "\n",
    "N += 1\n",
    "dx = torch.arange(N+0.0,device=device)*(dxmax-dxmin)/(N-1.0)+dxmin\n",
    "\n",
    "const = Constants().to(device)\n",
    "\n",
    "species = torch.as_tensor([[8,6,1,1]],dtype=torch.int64, device=device) \\\n",
    "               .expand(N,4)\n",
    "\n",
    "coordinates_op = torch.tensor([\n",
    "             [\n",
    "              [0.014497983896917479, 3.208059775069048e-05, -1.0697192017402962e-07],\n",
    "              [1.3364260303072648, -3.2628339194439124e-05, 8.51016890853131e-07],\n",
    "              [1.757659914731728, 1.03950803854101, -5.348699815983099e-07],\n",
    "              [1.7575581407994696, -1.039614529391432, 2.84735846426227e-06]\n",
    "             ],\n",
    "             ], device=device)\n",
    "#\n",
    "coordinates = coordinates_op.expand(N,4,3).clone()\n",
    "coordinates[...,at,:] += dx.unsqueeze(1)*dir.unsqueeze(0)\n",
    "coordinates.requires_grad_(True)\n",
    "\n",
    "#\"\"\"\n",
    "#f=open('setup.xyz','w')\n",
    "#for mol in range(N):\n",
    "#    f.write(\"%d\\n\\n\" % torch.sum(species[mol]>0))\n",
    "#    for atom in range(coordinates.shape[1]):\n",
    "#        if species[mol,atom]>0:\n",
    "#            f.write(\"%s %f %f %f\\n\" % (const.label[species[mol,atom].item()],\n",
    "#                                   coordinates[mol,atom,0],\n",
    "#                                   coordinates[mol,atom,1],\n",
    "#                                   coordinates[mol,atom,2]))\n",
    "#\n",
    "#f.close()\n",
    "#\"\"\"\n",
    "\n",
    "#np.save('test-R.npy', coordinates.detach().cpu().numpy())\n",
    "#np.save('test-Z.npy', species.detach().cpu().numpy())\n",
    "\n",
    "elements = [0]+sorted(set(species.reshape(-1).tolist()))\n",
    "seqm_parameters = {\n",
    "                   'method' : 'PM3',  # AM1, MNDO, PM#\n",
    "                   'scf_eps' : 1.0e-6,  # unit eV, change of electric energy, as nuclear energy doesnt' change during SCF\n",
    "                   'scf_converger' : [2,0.0], # converger used for scf loop\n",
    "                                         # [0, 0.1], [0, alpha] constant mixing, P = alpha*P + (1.0-alpha)*Pnew\n",
    "                                         # [1], adaptive mixing\n",
    "                                         # [2], adaptive mixing, then pulay\n",
    "                   'sp2' : [True, 1.0e-5],  # whether to use sp2 algorithm in scf loop,\n",
    "                                            #[True, eps] or [False], eps for SP2 conve criteria\n",
    "                   'elements' : elements, #[0,1,6,8],\n",
    "                   'learned' : [], # learned parameters name list, e.g ['U_ss']\n",
    "                   #'parameter_file_dir' : '../../seqm/params/', # file directory for other required parameters\n",
    "                   'pair_outer_cutoff' : 1.0e10, # consistent with the unit on coordinates\n",
    "                   }\n",
    "\n",
    "\n",
    "import time\n",
    "t0 = time.time()\n",
    "with torch.autograd.set_detect_anomaly(True):\n",
    "    eng = Energy(seqm_parameters).to(device)\n",
    "    Hf, Etot, Eelec, Enuc, Eiso, EnucAB, e, P, charge, notconverged = eng(const, coordinates, species, learned_parameters=dict(), all_terms=True)\n",
    "    L=Etot.sum()\n",
    "    L.backward()\n",
    "t1 = time.time()\n",
    "print('Time : %f sec' % (t1-t0))\n",
    "#print(Etot)\n",
    "#print(coordinates.grad)\n",
    "force = -coordinates.grad\n",
    "#f=open('log.dat', 'w')\n",
    "#f.write(\"#index, dx (Angstrom), energy (eV), force(eV/Angstrom)\\n\")\n",
    "#for i in range(N):\n",
    "#    f.write(\"%d %12.8e %12.8e %12.8e\\n\" % (i,dx[i],Etot[i],(force[i,at,:]*dir).sum() ))\n",
    "#f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.25\n",
      "-0.25\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Force')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAw3klEQVR4nO3dd3gU5dfG8e9JCBB6VZEWqiASQu+94w9EaSJFbBFRKQLSm6ICIkVBilIVG4qIiqBUpXeQ3nsHpQZIyHn/yMKLmLJCdifJns91zeXuzOzM/YDs2Zln5hlRVYwxxvgeP6cDGGOMcYYVAGOM8VFWAIwxxkdZATDGGB9lBcAYY3yUFQBjjPFRVgCMMcZHWQEwxkVEDopImIhcvmN62OlcxniKFQBj/qmhqqa5Yzru7gdFJJkngxkT36wAGBMLEUkhIqNE5LhrGiUiKVzLqonIURHpISIngSki4i8ivUVkn4hcEpH1IpLTtX4hEflNRM6LyC4Rae5o44zPswJgTOz6AOWAEKAYUAboe8fyh4BMQG4gFHgDaAk0ANIBzwNXRSQ18BvwBfCAa52PRaSIV1phTDTExgIyJoqIHASyABGuWUuAosDrqjrXtU5dYIKqBolINeBXIJ2qXnMt3wW8qao/3LXtFsBrqlr5jnkTgOOqOsiDzTImRnbO0ph/aqyqC269EZEw4NAdyw8Bd3YMn7n15e+SE9gXzXZzA2VF5O875iUDPrvvxMbcIysAxsTuOFFf3ttc73O55t1y9yH0ESAfsDWa+UtVtbYnQhpzL6wPwJjYfQn0FZGsIpIF6A98Hsv6nwJvi0gBiRIsIpmBn4CCItJGRAJcU2kRKeyFNhgTLSsAxsRuMLAO2AL8CWxwzYvJCOAbovoGLgKTgEBVvQTUAZ4m6gjiJDAUSOGx5MbEwTqBjTHGR9kRgDHG+CgrAMYY46OsABhjjI+yAmCMMT4qUd0HkCVLFg0KCnI6hjHGJCrr168/q6pZ757vWAFwDZA1naixVCKBiao6OrbPBAUFsW7dOm/EM8aYJENEDkU338kjgAigq6puEJG0wHoR+U1VtzuYyRhjfIZjfQCqekJVN7heXwJ2ANmdymOMMb4mQXQCi0gQUBxYHc2yUBFZJyLrzpw54/VsxhiTVDneCSwiaYDvgM6qevHu5ao6EZgIUKpUKbtt2RiHhIeHc/ToUa5duxb3ysYRKVOmJEeOHAQEBLi1vqMFQEQCiPryn6Gqs5zMYoyJ3dGjR0mbNi1BQUGIiNNxzF1UlXPnznH06FHy5Mnj1mccOwUkUf8HTQJ2qOoIp3IYY9xz7do1MmfObF/+CZSIkDlz5v90hOZkH0BFoA1QQ0Q2uaYGDuYxxsTBvvwTtv/69+PYKSBVXQZ45/+mX36BXbugaVPIkcMruzTGmIQuQVwF5HE//wxdukDOnFChAn+8/zpHdvzrgiNjTAL34YcfUrhwYVq1auV0lBhVq1Yt0dyw6hsFYMyYqCOAd95Br16h7bEx5PqmHBW6pGP0sKc4vnej0wmNMW74+OOPmTt3LjNmzPjH/IiICIcSJW6+UQAAChaE3r2RTZv5rc2vvOtXmzDC6Rz2PTk+L8Fbz+WBiRPh3DmnkxpjotG+fXv2799Po0aNGDlyJAMHDiQ0NJQ6derQtm1bDh06RM2aNQkODqZmzZocPnwYgHbt2vHKK69QvXp18ubNy9KlS3n++ecpXLgw7dq1i3Zfa9eupUKFChQrVowyZcpw6dIlrl27xnPPPUfRokUpXrw4ixcvBiAsLIynn36a4OBgWrRoQVhY2O3t/Prrr5QvX54SJUrQrFkzLl++7PE/p//C8fsAnJC/ZG16laxNL2DX6rl8/cv7VN63D6a+zM7+HXijZSZaFnyKxi0GkDZTNqfjGpPwdO4MmzbF7zZDQmDUqBgXjx8/nnnz5rF48WKyZMnCwIEDWb9+PcuWLSMwMJCGDRvStm1bnn32WSZPnkzHjh2ZPXs2AH/99ReLFi1izpw5NGzYkOXLl/Ppp59SunRpNm3aREhIyO393LhxgxYtWvD1119TunRpLl68SGBgIKNHRw1V9ueff7Jz507q1KnD7t27GTduHKlSpWLLli1s2bKFEiVKAHD27FkGDx7MggULSJ06NUOHDmXEiBH0798/fv/c7oPvHAHE4JGyDeg/cDHVlx6CjRs58lwTtiU7T9vTE3hwxMO07Jqbn78YSPi1q05HNcbcpVGjRgQGBgKwcuVKnnnmGQDatGnDsmXLbq/XsGFDRISiRYvy4IMPUrRoUfz8/ChSpAgHDx78xzZ37dpFtmzZKF26NADp0qUjWbJkLFu2jDZt2gBQqFAhcufOze7du/n9999p3bo1AMHBwQQHBwOwatUqtm/fTsWKFQkJCWHatGkcOhTtmGyO8ckjgGiJQEgItUO+5sDNCFbMncCMPz7mm4AdfL9jECfzjyFDk1acb9GIjOWqI34+XzuNL4vll7o3pU6dOsZld14SmSJFCgD8/Pxuv771/u7+A1WN9nLK2J6fHtP6tWvX5ssvv4y5AQ6zb7Fo+Pkno1LDVxk3bBsnBlzkj0eHk6F8dRg/nvpTalGkWyBD3qnPsd2Jo6ffGF9QoUIFvvrqKwBmzJhBpUqV7mk7hQoV4vjx46xduxaAS5cuERERQZUqVW53Pu/evZvDhw/zyCOP/GP+1q1b2bJlCwDlypVj+fLl7N27F4CrV6+ye/fu+2pjfLMCEIfkgWko3bIrzJyJnjjBi0VakykyBb0i5pFrRmnqd8nKoin94fp1p6Ma49M+/PBDpkyZQnBwMJ999tntc/b/VfLkyfn66695/fXXKVasGLVr1+batWt06NCBmzdvUrRoUVq0aMHUqVNJkSIFr7zyCpcvXyY4OJhhw4ZRpkwZALJmzcrUqVNp2bIlwcHBlCtXjp07d8Znk++bxHZYk9CUKlVKE8r1tXs3LGDa7IFMC1tFv0U3eelgJi60acahprUJrtTE6XjGxLsdO3ZQuHBhp2OYOET39yQi61W11N3r2hHAPcpfohZvv7WMA+9e5dl3foJatfhszacUW9iUMl3S8Mnotlw+f9LpmMYYEyMrAPfJPyA5yes9Dl9/TctvdjAq8EmuSgShf3/Gw8Oz8UqPIkRsSBhHLcYYcycrAPEoc44CdHpzFn8Ov8ry0uN56kY+Dp7YSbKSpaOGoJjYl2uX/3Y6pjHGAFYAPEL8/KjQ4GWmDt/L3JGnYcQITl8+Rc2j75BzcCZ696vA4e0rnY5pjPFxVgA8TDJnhi5dyLppD3ODh1Lp+kMM9VtJnq8r0OSN7OyeNwMSUUe8MSbpsALgJeLnR62mb/L9yOPsb76M7jfL8nvACZK3bA0lS3Ji8ofcCEtY44QYY5I2KwAOyF2kIkMGr+J4378IGjoBrl/nxYWdCBqQnvcG1+X88X1ORzTGRGPOnDkMGTLknj4bFBTE2bNnY11n6tSpvPbaa0DU2EfTp08HYOfOnYSEhFC8eHH27dsXb8NiWwFwUEDa9BAaClu30rHhYB67kYHeN38l19j8dO5VnEPbljsd0RjjEhERQaNGjejZs6dX9te+fXvatm0LwOzZs3niiSfYuHEj+fLli3FY7P/KCkBCIELdp/vw66hzbK4xk6eu52VswCY+e7UytGkDW7c6ndAYxx08eJDChQvz0ksvUaRIEerUqXN76OU7H8Jy9uxZgoKCgKhf1I0bN6Zhw4bkyZOHMWPGMGLECIoXL065cuU4f/48APv27aNevXqULFmSypUr375jt127drzxxhtUr16dHj16/OMX+qlTp3jyyScpVqwYxYoVY8WKFQA0btyYkiVLUqRIESZOnBhnu6ZMmULBggWpWrUqy5f//4++gQMHMnz4cObOncuoUaP49NNPqV69+r+Gxb4fNhhcAhNcuSnTKzflnR2rSXt5GkyczrcbPufzOg/Su967lKn7vNMRjQGg2tRq/5rXvEhzOpTuwNXwqzSY8e9HfLcLaUe7kHacvXqWpt80/ceyJe2WxLnPPXv28OWXX/LJJ5/QvHlzvvvuu9sjccZk69atbNy4kWvXrpE/f36GDh3Kxo0b6dKlC9OnT6dz586EhoYyfvx4ChQowOrVq+nQoQOLFi0Cosb9WbBgAf7+/kydOvX2djt27EjVqlX5/vvvuXnz5u2x/idPnkymTJkICwujdOnSNGnShMyZM0eb7cSJEwwYMID169eTPn16qlevTvHixf+xToMGDWjfvj1p0qShW7duAP8YFvt+2BFAApWzcFkyjPgYDh3iQvNG/JHiNGVXvUCdzpn5Y85HTsczxhF58uS5PXZ/yZIl/zWUc3SqV69O2rRpyZo1K+nTp6dhw4YAFC1alIMHD3L58mVWrFhBs2bNCAkJ4eWXX+bEiRO3P9+sWTP8/f3/td1FixbxyiuvAODv70/69OmBqDGJihUrRrly5Thy5Ah79uyJMdvq1aupVq0aWbNmJXny5LRo0cLdP4p4YUcACV3mzLww4AdanD/JuPEvMDzFL1TZ2JF2c99hytNfQbVqTic0Piq2X+ypAlLFujxLqixu/eK/251DOfv7+98+BZQsWTIiIyMBuHbtWoyfuXM46FtDQUdGRpIhQwY2xfCAm9iGnL7bkiVLWLBgAStXriRVqlRUq1btX3nuFt1Q0t5iRwCJRJpMD9G9988c6H2aUYFP0mDrdahenas1q7Bszhin4xnjqKCgINavXw/At99++58+my5dOvLkycPMmTOBqHH8N2/eHOfnatasybhx4wC4efMmFy9e5MKFC2TMmJFUqVKxc+dOVq1aFes2ypYty5IlSzh37hzh4eG3M3iLFYBEJlX6LHR6cxbNfjsOo0cz2W8zlTe+Tt3OWVgzf7LT8YxxRLdu3Rg3bhwVKlSI81LL6MyYMYNJkyZRrFgxihQpwg8//BDnZ0aPHs3ixYspWrQoJUuWZNu2bdSrV4+IiAiCg4Pp168f5cqVi3Ub2bJlY+DAgZQvX55atWrdfpykt9hw0Inc1YvnGPfxcwz5+yfOBipP/P0g7zQdR5GKTzodzSQxNhx04mDDQfuQVOky07XnHPZ3P8rbfjVZHHiK9pOegrZt4cABp+MZYxIwKwBJRNrMD9O33wL2v7qHSQ+8CDNncqrkI/TsVZrzx/Y6Hc8YkwBZAUhiMmfPT8Ehn8DevfzWpgLDUqwj35gCfPDO44RfueR0PJPIJaZTxr7ov/79WAFIqrJnp/XoJWyuM4vy17LSLWIupXpnZvWUwTb6qLknKVOm5Ny5c1YEEihV5dy5c6RMmdLtz1gnsI+Y/c1bvLZhMJX2hfPV0XIwahSULet0LJOIhIeHc/To0TivazfOSZkyJTly5CAgIOAf82PqBLYC4EMuXf2b619MJ0u/99gaeZLFzUvTvsdMAnLkdjqaMcaDEuRVQCJST0R2icheEfHOEHs+LG2qDGR5sSPs3s2M0HJ0zLKWYsPyMOu9tuj1607HM8Z4mWMFQET8gbFAfeBRoKWIPOpUHp+SNi3vvrWC76uOQ1OnpsmNzyjTPT0Lv7m3cc6NMYmTk2MBlQH2qup+ABH5CngC2B7fOxr04za2H78Y35tNAopRNGg+D5zaz6HAY/RaouRZ/inkyw93jJ9ijHHeow+nY0DDIvG6TScLQHbgyB3vjwL/6pUUkVAgFCBXrlzeSeZDBOGhB/PxYGQQevQoHDrMX1dXc/mBjOTI8Sji9+9REI0xSYOTBSC6IfD+1SOtqhOBiRDVCXwvO4rvqpmkHTjA68NrMCb5QfLuT86wigN56omejo5YaIzxDCc7gY8COe94nwM47lAWc0uePHw09gDz8vQn8EYkTTf3plrf7Ow4sMbpZMaYeOZkAVgLFBCRPCKSHHgamONgHnOHum0Hsemt04y/UoOtESdY9mIdcGOERGNM4uFYAVDVCOA1YD6wA/hGVbc5lcf8W7L0GXl52EJ2N/qVF8/mgsaNmfV8eTZtXeB0NGNMPLAbwYx7wsOJHDaUx072Z18GZcgDz9Dp1en4WSexMQlegrwRzCQiAQH49enL7y8up975jLxx/gsadH2IY7vWOp3MGHOPrACY/yRLsfLMHn2acama8Ueqszw2tQwnJ422AeaMSYSsAJj/TJIlo333b9jcbAF9j+ThoRc7Q6NGXDtiD6AxJjGxAmDuWf6QmnSdvhdGjWLjn7+S56N8zPy0i9OxjDFusgJg7o+fH3TqROqZP5AzPBXNj43imTdy89fx/U4nM8bEwQqAiRcFS9djxdCzvCU1mJnmMMEjC7L4uw+cjmWMiYUVABNvkiVPSb/+C1lZZRqpIv1YPLYbdO0KNtS0MQmSFQAT70rVasuGfsfoV/hlGDGC5fUfY8eqH52OZYy5ixUA4xGpM2QlYOx4dM4cOhQ5QMmfGjF+xDNoZKTT0YwxLlYAjEdJw4bM67yeSlcy88qlL2nRNRcXTx50OpYxBisAxguy5SvGvOEnGZL8cWalO0bpoQU4tdDG/TPGaXEWABEpJSJdROR9EXlLRJqLSCZvhDNJh59/Mnr0+olFFSZQ7UwqHqjTGN55B+yUkDGOibEAiEg7EdkA9AICgV3AaaAS8JuITBMRe0SX+U+q1A1lwrijyNMtOTi8L21fz8G5QzucjmWMT4rtiWCpgYqqGhbdQhEJAQoAhz2QyyRladPC55+z+qM0fHV2IgvHPsZn5d+nxpNvOJ3MGJ8S4xGAqo6N6cvftXyTqi70TCyT5InQouMEVtf6irQ3k1Frc1f6vF2N8HC7Z8AYb3GnDyCPiIwQkVkiMufW5I1wJukrXqUF6/sc4vkLeXk3cinD2heFs2edjmWMT3DnKqDZwEHgI+CDOyZj4kXqTA/x6Yi9fJ/+ZTp9fRBKluTmmlVOxzImyXOnAFxT1Q9VdbGqLr01eTyZ8S0iNO48njRLVnApQCk7uQJTRz5rzxkwxoPcKQCjRWSAiJQXkRK3Jo8nM76pVCkili4mQ6pMPHdxOp26PsqNyxecTmVMkuROASgKvAQM4f9P/wz3ZCjj2zJmz8cvQ4/RiXJ8mH4nVfs8zJGty52OZUyS404BeBLIq6pVVbW6a6rh6WDGtwUEpGDUgJV8k68n21Jf5aUPqsH8+U7HMiZJcacAbAYyeDiHMdFq1vo91jX7jfE78kGDBlwd8jaRkTedjmVMkuBOAXgQ2Cki8+0yUOOEgsVrEbRwPdq0Cc9s7k/DN3Nw/ozdf2jM/YrtTuBbBng8hTFxSZ0avvyKukPC6XRtNmWHFeDHVj9TKKSW08mMSbTiPAJwXfK5E0jrmnbYZaDGCeLnxyu9v2dJsRFc8Aun3Dd1+O27YU7HMibRcudO4ObAGqAZ0BxYLSJNPR3MmJhUaNKFtS0XkyssOe1W9CBs9HC7X8CYe+DOKaA+QGlVPQ0gIlmBBcC3ngxmTGxyh1RleY59HHy1NYEjuhO5bReRH44mWcpUTkczJtFwpxPY79aXv8s5Nz9njEelzZKdol8uhN696b/vU+q/mZ2/juxxOpYxiYY7X+TzXFcAtRORdsDPwFzPxjLGTX5+8M475GsWytIMf1Nu5KPsWfmz06mMSRRiLQAiIsCHwAQgGCgGTFTVHl7IZozbnms/gYUVJ3Au+U3K/vA/Fn/5rtORjEnwROPoPBOR9apa0kt5YlWqVCldt26d0zFMArZ/xwoaTqrJwRTXOJhxIFm79gcRp2MZ4yjX93ipu+e7cwpolYiUjucw74vIThHZIiLfi0iG+Ny+8V15C1dgRZ/9zDpeiazdB0JoKJHXrzkdy5gEyZ0CUB1YKSL7XF/Yf4rIlvvc72/AY6oaDOwm6rnDxsSL9BmzUXfSUujbl69XfkrNntk4c2y307GMSXDcuQy0fnzvVFV/vePtKsDuKzDxy88P3n6byIfPserYOMqOeoy5T/9IoZJ1nU5mTILhzhHAYFU9dOcEDI7HDM8Dv8S0UERCRWSdiKw7c+ZMPO7W+IKWr3zM0vLjueIXQflv67Pkx4+cjmRMguFOAShy5xsR8Qfi7BQWkQUisjWa6Yk71ukDRAAzYtqOqk5U1VKqWipr1qxuxDXmn8o8/jKrWi4k2/UA6qztyO4p9jgLYyCWU0Ai0gvoDQSKyMVbs4EbwMS4NqyqsY7SJSLPAv8DampclyIZc5/yhFRn+QM7+LZrfQq+3R2OXYM+fewKIePTYjwCUNX3VDUt8L6qpnNNaVU1s6reV6etiNQDegCNVPXq/WzLGHdlfDgvL037E9q0Yf3H/ejQrTA3wi47HcsYx7hzCugnEUkNICKtRWSEiOS+z/2OIWpk0d9EZJOIjL/P7RnjnuTJYdo0fn+5LuPS7aJBr1z8feqQ06mMcYQ7BWAccFVEigFvAoeA6fezU1XNr6o5VTXENbW/n+0Z85+I0GXAPKZmeYnf0/1FxaGPcHCbPXPY+B53CkCE6xz9E8BoVR1N1K93YxK1Z1+dyPyQ4RxPcZ2y06uw948fnI5kjFe5UwAuuTqEWwM/u64CCvBsLGO8o/pTXVnZaA6NDwcS1OAZmDfP6UjGeI07BaAFcB14QVVPAtmB9z2ayhgvKlS+IRM+2E2y/AU51eJxZnwU6nQkY7zCnUdCnlTVEar6h+v9YVW9rz4AYxKchx+G339n2NM5aX3+EwYMrIpGRjqdyhiPcueRkJdE5OJd0xHXIG55vRHSGK9Im5Yho7fT7nJ+3pLfCe1RmIjrYU6nMsZj3DkFNALoTtSpnxxAN+AT4CtgsueiGeN9ASlTMXnoLvpEVuLTNLt5qkcQV/86HfcHjUmE3CkA9VR1gqpeUtWLqjoRaKCqXwMZPZzPGK8TPz8GD/qDselasi/8NGEN6sCpU07HMibeuVMAIkWkuYj4uabmdyyzIRxMktWhyxdsqDuLzJt3c71yBQ5v/t3pSMbEK3cKQCugDXAaOOV63VpEAoHXPJjNGMelaPQkLF5Ml6LHKP1FNdb9Ns3pSMbEG3euAtqvqg1VNYuqZnW93quqYaq6zBshjXFU2bJ06jWbVDf9qbqkHb98+ZbTiYyJF+5cBZRVRHqLyEQRmXxr8kY4YxKKR0rVY2WH9TxyNZCGOwcwbexLTkcy5r65cwroByA9sAD4+Y7JGJ/yUN5glvTeTfULmeh85FPOvdcPbCRzk4i580jIVKraw+NJjEkE0mXNwc/vHWLXK83JPHQwevICOmIEfv7u/FMyJmFxdzjoBh5PYkwikTwwDUUn/wRdujB83Ue06p6P61cuxv1BYxIYdwpAJ6KKQJjrLuBLdzwhzBjf5OcHH3yAPP44X6U/zON9grh45qjTqYz5T9y5CiitqvqpauAdTwVL541wxiRoInTr/RPTsrzEkvR/Ue3dgpw88KfTqYxxmztHALeJSD4R6SMiWz0VyJjEpu2rE/mx0CB2pQqj8piSXNu93elIxrjFnctAs4lIFxFZA2wjquO4pceTGZOI1G/Zn8VVJtFrTXJSVqkBGzc6HcmYOMVYAETkJRFZBCwFMgMvAidUdZCq2nGuMXcpU/d5nv9kHaRIwfy2FVny/UinIxkTq9iOAMYC/sAzqtpXVbdgY/8YE7tChYhc9gf9akC9DW/wwxS7gtokXLEVgIeJGvJ5hIjsEpG3sUdBGhMnv5y5+OXNLRS7koYmB4cxeWRbpyMZE60YC4CqnlXVcapaBagJXABOi8gOEXnXawmNSYQyZ8/Pwv77qHkpCy9c/Ix3B9WyJ4yZBCe2PoBst16r6lFVHa6qJYHGRD0j2BgTizQZHuDHdw/S6nJeDqxfCB1fh5s3nY5lzG2xnQKaLCKrRGSIiFQTkWQAqrpLVQd5KZ8xiVrylKmZPnQ34wt1RcZ+zOG2T9hdwybBiO0UUH2gGrAEeBJYJSKzRCRURHJ5J54xiZ+fnz/+w4YTNuxdqmX9mQZ9g7h49pjTsYyJ/T4AVb2mqvNUtZOqlgK6EnUfwBjXfQHGGDcFdu/FWwVC+T3tX1R9twCnDm1zOpLxcbH1AYwRkYp3zlPVA6r6sao2Aip5PJ0xSUzrVyfw4yMD2R0YRqWPinNwqz1TyTgntiOAPcBwETkoIkNFJOTOhap6w6PJjEmi6rUawIIK4zkbEEG392vDVhtZxTgjtj6A0apaHqgKnAemuC4B7S8iBb2W0JgkqPzjL7Os4Sw++T0DVK4My5c7Hcn4IHdGAz2kqkNVtTjwDFEdwjs8nsyYJK5IhcZkXLKKaw9lodG4Kvz6hT1r2HiXO4PBBYhIQxGZAfwC7AaaeDyZMb4gd24uzvuBQw+m4H87B/Dl2PZOJzI+JLZO4Nquh78fBUKBuUA+VW2hqrPjY+ci0k1EVESyxMf2jEmMHsj9KEt77qL85Qw8c3YCH777hNORjI+I7QigN7ASKKyqDVV1hqpeia8di0hOoDZwOL62aUxilSFrTua/fZAnL2anU/gchvWuAjZ0hPGw2DqBq6vqJ6p6XkQqichzACKSVUTyxMO+RwJvYiOMGgNAytTpmTn0AF2vhvD4p39Au3YQHu50LJOEudMHMADoAfRyzQoAPr+fnYpII+CYqm52Y91QEVknIuvOnDlzP7s1JsHzTxbA8CEbKNLxbfSzz5j4Ugmu/HXa6VgmiRLV2H+Ai8gmoDiwwXUlECKyRVWD4/jcAuChaBb1Ier0Uh1VvSAiB4FSqno2rrClSpXSdevWxbWaMUnChrF9KX36HcpcTMPP3TaSKXt+pyOZREpE1rtGc/gHd54JfEOjqoS6NpTanR2qai1VfezuCdgP5AE2u778cwAbRCS6YmGMzyrx6mBm5nmTDWkuU+mDIhzZvsrpSCaJcacAfCMiE4AMIvISsAD45F53qKp/quoDqhqkqkFEXWVUQlVP3us2jUmqnmo3lPmlRnEs5Q3KT6nIthU/OB3JJCHu3Ag2HPgW+A54BOivqh95OpgxJkq1Jzrxe/1vEGDf663srmETb2LsAxAR0Tg6CNxZJz5ZH4DxZWF7dxL4+BNw+DBHPhtDzqYvOB3JJBL30gewWERev3vsfxFJLiI1RGQa8Gx8BzXGRC8wfyFYtowl1YLIv/lFxo1q5XQkk8jFVgDqATeBL0XkuIhsF5EDRI0S2hIYqapTvZDRGHNL1qyU/mIpdS5lpcOFL+gzqIo9a9jcszgvA4Wo8YCALECYqv7t6VAxsVNAxkSJuB5Gh97F+CTdHtpdLcgng7eQLCCF07FMAnU/l4GiquGqesLJL39jzP9LliKQCe/vZGB4Raam2s3M0EoQFuZ0LJPIuFUAjDEJj/j5MWDwMpam7cjT09ZB7dpw/rzTsUwiYgXAmESuyhujkW9msmP/Gqr3zcmxHfa4buMetwqAiOQWkVqu14EiktazsYwx/0nTppz6eBjrM1ylwqQK7Fj2vdOJTCLgzmBwLxF1I9gE16wcwGwPZjLG3INqjTuztME3XPdXKv30FCtmjXY6kkng3DkCeBWoCFwEUNU9wAOeDGWMuTfFKzVjxbO/kzkiOTU2dmb5J/2cjmQSMHcKwHVVvXHrjYgkw8bwNybByvtoRVZ03U77kzko2WEwvPceeO+GfZOIuFMAlopIbyBQRGoDM4EfPRvLGHM/smTLx6ix+0jZohV/vdWbQd1LE3HjmtOxTALjTgHoCZwB/gReJurZwH09GcoYEw+SJ4fp05ndtQED066nUY9cXD5vg+6a/+dOAQgEJqtqM1VtCkx2zTPGJHR+fjw3+GcmpGvFr+nOUPXtvJzYt8npVCaBcKcALOSfX/iBRD0TwBiTSIR2+ZwfHxnIrlRhlB1fiu0r5zgdySQA7hSAlKp6+dYb1+tUnotkjPGE+s8M4I+an5P9sh8Zm7aGRYucjmQc5k4BuCIiJW69EZGSgA06YkwiVLxGK1b02E22jLmIqFeH2R93dDqScZA7BaATMFNE/hCRP4Cvgdc8G8sY4ykSFATLlzOleUGePPMRnfqX4WZEuNOxjANiLQAi4g9UBgoBrwAdgMKqut4L2YwxnpI+Pc9P2kCXS4/xof9aGvcM4vLFs06nMl4WawFQ1ZvAE67hoLe6HuhuPxWMSQL8U6RkxPtbGBvQmLmpj1NlUBDHD2xxOpbxIndOAS0XkTEiUllEStyaPJ7MGON5InTo/T0/5enNcb8rHG9WD7ZvdzqV8ZI4nwgmIoujma2qWsMzkWJmTwQzxnOurvqDVE82h6tX2TNtJAUaP+90JBNP7vmJYKpaPZrJ61/+xhjPSlWuMqxZwzeVMlJ44wtMGP6005GMh7kzHHR6ERkhIutc0wcikt4b4YwxXpYzJ/U/W0ndSw/Q/srX9HizBJE3rjudyniIO30Ak4FLQHPXdBGY4slQxhjnpM2UjR+GHqZDeHGGpd5I8645uXLqiNOxjAe4UwDyqeoAVd3vmgYBeT0dzBjjnGQBKRjz9npGZGzJ95nOsLhFWdixw+lYJp65UwDCRKTSrTciUhG7E9iYJE9E6NLxC3ZV/pr/7bgJ5cpxfs7XTscy8cidAtAeGCsiB0XkIDCGqGGhjTE+IH+t5rB2LatKPkjuVU8z7t0n0chIp2OZeBBjARCRXACqullViwHBQLCqFldVu1vEGF+SKxePfLuYKuHZ6BA+m+e65SfswjmnU5n7FNsRwOxbL0TkO1W9qKoXPR/JGJMQZcyUnR+HHmGAXw2mpT9AxYE5ObhtudOxzH2IrQDIHa+t09cYg5+fPwP7LeSn/APYnzKMz7vVgcXR3StqEoPYCoDG8DpeiMjrIrJLRLaJyLD43r4xxnMebzWQLc0X0+twLqhdm6MfDLB+gUQotgJQTEQuisglINj1+qKIXBKR+zoVJCLVgSeI6lMoAgy/n+0ZY7wvV/Fq+K9aw+mn6lLy5Fu07J6Hy3+dcjqW+Q9iLACq6q+q6VQ1raomc72+9T7dfe73FWCIql537ev0fW7PGOOEtGnJ+uUcuqSvy8w0hyn3dm72rP/N6VTGTe5cBuoJBYHKIrJaRJaKSGmHchhj7pP4+9Oz7zzmPzaEk8lvUOq7Osye1svpWMYNHisAIrJARLZGMz0BJAMyAuWA7sA3IiIxbCf01jhEZ86c8VRcY8x9qtWsB+tb/84jYamY/ssQtMebEBHhdCwTiziHg/bITkXmEXUKaInr/T6gnKrG+g1vw0Ebk/Bdv3KRG927kHbcZA7VKUvA2HE8nL+407F82j0PB+0hs4EaACJSEEgO2PPojEkCUqROR9qPJ8H06TybYy3FPynJwu/edzqWiYZTBWAykFdEtgJfAc+qE4cixhjPadOGj1/8nizhAdT+800GDarOzfAbTqcyd3CkAKjqDVVtraqPqWoJVV3kRA5jjGc9Wr4Ra/odps3lvAxkCXW6P8SZ/VudjmVcnDoCMMb4iNQZH2TqsD1MzvQcl8L+JnXlGvDrr07HMlgBMMZ4gfj58dzrk1n1+mZSZXyAyw3rMrRPNW6EXXY6mk+zAmCM8Rq/x4rCmjXM7lCdnsmXUrH3g+zduNDpWD7LCoAxxrtSpaL1yEV8l6s7+1KGUfzbWkwd86KNJeQAKwDGGEc89dwwNrddSckr6Xnu3CTefq0onD/vdCyfYgXAGOOYnIXLsvD907zvX582M3dBcDARv813OpbPsAJgjHGUf0ByuvWdS555q9E0qXlqSj069izG1Qt2b6inWQEwxiQMJUsSsXYNeXIX46PALZQclJ018yc7nSpJswJgjEkwAtKmZ/R7m/ityFAu+9+kwooX6NOvItev2NNoPcEKgDEmwanV9E22dtvHs1cLMDlsBZerlIX1652OleRYATDGJEjpH8zNpPd3s7XyF2Q+cYHw8mUY3bcW1y7/7XS0JMMKgDEmQcv8REvYvp1fX6pJ54CFFO//ICt+Hu90rCTBCoAxJuHLkIHHx/7K/MLvctX/JpXWvsKrPR7j4ukjTidL1KwAGGMSjTrNe7Gt11E6XS/OuMBtNO9dAGbNAhtN/p5YATDGJCppMj3EyPc2sKriFN7ZmwuaNOHvJ+tzaNtyp6MlOlYAjDGJUpna7Si5YDsMH86AyIU8+kUlhg6ubyOM/gdWAIwxiVeyZNC1K12H/E7tsIfoeXMexfpkZuHMYU4nSxSsABhjEr1cj5Zn9ogT/FxwEOESSa3tPfjg5WA4eNDpaAmaFQBjTJLRoGV/tg46w9t+NWn80x4oVIhT/bpw9e8zTkdLkKwAGGOSlJRpMtC33wLyrd4DTZrw4r5RFHznIT4bE0rkzQin4yUoVgCMMUlTjhwwYwY9Wn1MtohA2p77hDLd0rH0+1FOJ0swrAAYY5K0So+/wurhf/P5Qx04lewG1bZ0YdKzRWHLFqejOc4KgDEmyfPzT0arl8eyu/8ZhgU8zlPzD0NICGtDH2fnpgVOx3OMFQBjjM8ITJuR7r1/IuP2A1GXj0b+wqOza9OsV342b/W9QmAFwBjjezJlgvff57uem+h1pTi/so+Q72rzVJ/8/LltsdPpvMYKgDHGZ2XNH8w772/g0LMbGfB3CAsj97Hyxbrwxhtw/LjT8TzOCoAxxudlKBTCwJEbOdh6He0KNocPP2TCU7l4rmchDm5d5nQ8j7ECYIwxLhmLlCT5tM9h927OVyzBlwG7yD+zMq275WHLsu+cjhfvrAAYY8zd8ual1wdr2Nt6DR3DSzA7xUGKLWxK79B8sHBhkhl+2gqAMcbEIMcjpRnx7noOv7qXt6UmVTecg1q1OFnuMaaNfSnRP57SCoAxxsQh08P56Nt/AXWXnYRJk/jy4XO0O/spOQdnok//ihzdtdbpiPfEkQIgIiEiskpENonIOhEp40QOY4z5T1KmhOefp/N3x1nw2DAqXn+Q9/xWEPRFGZp1yU7kvF8gMtLplG5z6ghgGDBIVUOA/q73xhiTKIifHzWbdGf2yBPse2op3SPLkfHYefzqN4D8+Zn+VpNEcVSQzKH9KpDO9To9kPQvuDXGJEl5gqvwXvBKuH4dnvqeE1M/4rmbs+CLWdS9mJV2RVrR6OkBpEyTwemo/yLqQG+2iBQG5gNC1FFIBVU9FMO6oUAoQK5cuUoeOhTtasYYk2Ds27SIKbP6My1sFUfT3CT9NZh1qQE1mr0JlSuDn3dPvojIelUt9a/5nioAIrIAeCiaRX2AmsBSVf1ORJoDoapaK65tlipVStetWxfPSY0xxjNuht9g8ZxRTF85gQ+mniDruTC+r5iZrVUe4enaXShQrQmIeDyH1wtAHGEuABlUVUVEgAuqmi6uz1kBMMYkWleuwI8/8sbvfRj54H4AQs4F0CxNWRpVDeWxWs+Av79Hdh1TAXCqE/g4UNX1ugawx6EcxhjjHalTw9NPM+LjfRx+bgsjM7YkMEVq+qRYRujMtpA9O4SGsuarD7h28bxXIjl1BFAJGE1UJ/Q1oIOqro/rc3YEYIxJao4f28np32YT8stGLv32M1levUKySBh48hG6F+sAdetCwYL3daooQZ0CuldWAIwxSVn41css/PkjftzwFbVXnabxkpNRC3LnhilToHr1e9puTAXAqctAjTHG3CUgVRrqNetFvWa9ombs3w/z50dNOXLE+/7sCMAYY5K4hNYJbIwxxmFWAIwxxkdZATDGGB9lBcAYY3yUFQBjjPFRVgCMMcZHWQEwxhgfZQXAGGN8VKK6EUxEzgCJ8YEAWYCzTofwIl9rL1ibfUVibXNuVc1698xEVQASKxFZF91deEmVr7UXrM2+Iqm12U4BGWOMj7ICYIwxPsoKgHdMdDqAl/lae8Ha7CuSVJutD8AYY3yUHQEYY4yPsgJgjDE+ygqAB4hIJhH5TUT2uP6bMZZ1/UVko4j85M2M8cmd9opIThFZLCI7RGSbiHRyIuv9EpF6IrJLRPaKSM9olouIfOhavkVESjiRMz650eZWrrZuEZEVIlLMiZzxKa4237FeaRG5KSJNvZkvvlgB8IyewEJVLQAsdL2PSSdgh1dSeY477Y0AuqpqYaAc8KqIPOrFjPdNRPyBsUB94FGgZTRtqA8UcE2hwDivhoxnbrb5AFBVVYOBt0nkHaVutvnWekOB+d5NGH+sAHjGE8A01+tpQOPoVhKRHMDjwKfeieUxcbZXVU+o6gbX60tEFb3s3goYT8oAe1V1v6reAL4iqu13egKYrlFWARlEJJu3g8ajONusqitU9S/X21VA/D+81rvc+XsGeB34DjjtzXDxyQqAZzyoqicg6osPeCCG9UYBbwKRXsrlKe62FwARCQKKA6s9Hy1eZQeO3PH+KP8uYu6sk5j81/a8APzi0USeF2ebRSQ78CQw3ou54l0ypwMkViKyAHgomkV93Pz8/4DTqrpeRKrFYzSPuN/23rGdNET9auqsqhfjI5sXSTTz7r6O2p11EhO32yMi1YkqAJU8msjz3GnzKKCHqt4UiW71xMEKwD1S1VoxLRORUyKSTVVPuA7/oztErAg0EpEGQEognYh8rqqtPRT5vsRDexGRAKK+/Geo6iwPRfWko0DOO97nAI7fwzqJiVvtEZFgok5l1lfVc17K5inutLkU8JXryz8L0EBEIlR1tlcSxhM7BeQZc4BnXa+fBX64ewVV7aWqOVQ1CHgaWJRQv/zdEGd7JepfyiRgh6qO8GK2+LQWKCAieUQkOVF/b3PuWmcO0NZ1NVA54MKt02OJVJxtFpFcwCygjarudiBjfIuzzaqaR1WDXP9+vwU6JLYvf7AC4ClDgNoisgeo7XqPiDwsInMdTeYZ7rS3ItAGqCEim1xTA2fi3htVjQBeI+qqjx3AN6q6TUTai0h712pzgf3AXuAToIMjYeOJm23uD2QGPnb9va5zKG68cLPNSYINBWGMMT7KjgCMMcZHWQEwxhgfZQXAGGN8lBUAY4zxUVYAjDHGR1kBMMYYH2UFwBhjfNT/AY3vC0rmQLzmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d = np.array([np.arange(0,101,1),dx.detach().cpu().numpy(),Etot.detach().cpu().numpy(),torch.sum((force[:,at,:]*dir),axis=1).detach().cpu().numpy()]).T\n",
    "#d=np.loadtxt(fn)\n",
    "##index, distance, energy, force\n",
    "\n",
    "t=np.argmin(d[:,2])\n",
    "print(d[t,1])\n",
    "\n",
    "\n",
    "t=np.argmin(np.abs(d[:,3]))\n",
    "print(d[t,1])\n",
    "\n",
    "#-0.122\n",
    "\n",
    "dx = d[1,1]-d[0,1]\n",
    "\n",
    "x1 = d[1:-1,1]\n",
    "f1 = -(d[2:,2]-d[:-2,2])/dx/2.0\n",
    "plt.figure(0)\n",
    "plt.plot(d[:,1],d[:,3],'r',label='from code')\n",
    "plt.plot(x1,f1,'g--', label='numerical diff')\n",
    "plt.plot(d[:,1], np.zeros_like(d[:,1]))\n",
    "plt.legend()\n",
    "plt.ylabel('Force (eV/Angstrom)')\n",
    "plt.title(\"Force\")\n",
    "#plt.savefig(\"force.jpg\")\n",
    "#plt.xlim([0.0,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Force difference')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEICAYAAACEdClSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABOVklEQVR4nO2deZwkVZXvvyf32tfed6BZpWmhWdxwARRQX7uMCjKCywzyRtTZxWUWZ/Q9fe/pqOOCiCg4KoM6aqvtIKCIgiyNNM3a0N303l1dXfuWWZWZ5/0REVlZWblVVWZWVdb5fj75yYyIeyNudFfmid+5554jqophGIZhVALfbA/AMAzDWDiY0TEMwzAqhhkdwzAMo2KY0TEMwzAqhhkdwzAMo2KY0TEMwzAqhhkdw5jDiMg/i8h/uJ9Xi8igiPjd7SUicp+IDIjI58ThWyLSIyIPz+7IDSM7gdkegGHMJiKyF1gCJNJ2n6yqh2dnRLlR1f1Afdqua4HjQKOqqoi8ArgEWKmqQ7MxRsMohCkdw4A3qmp92mtKBkdEZuvhbQ3wtI6v8F4D7J2OwZnFezAWGGZ0DCMLIhIWkS+IyGH39QURCbvHXiUiB0XkIyJyFPiWiPhF5GMistt1dz0qIqvc9qeKyF0i0i0iO0Xk7Xmuu05Efuue4y6gPe3YWhFREQmIyLeBa4C/d11u7wduBl7ibn/S7fMGEdkuIr0i8oCIbEg73173HnYAQ+55L3Db9YrI4yLyqrT294rIv4rI/e74fiUi6eN7eVrfAyLy7rR/y/8nIvtFpENEbhSRmhn/JxnzE1W1l70W7AvYC1ycZf+/AA8Ci4FFwAPAv7rHXgXEgc8CYaAG+DvgCeAUQICzgDagDjgAvAfHnX02jkvsjBzj+QPwefe8FwIDwH+4x9YCCgTc7W8Dn0rr+27g92nbZwPHgPMBP46R2guE0+59O7DKvYcVQBdwOc4D6SXu9iK3/b3AbuBkt/29wGfcY6vdsV4JBN173+ge+wKwBWgFGoCfAf97tv/v7TU7L1M6hgE/cZ/Oe0XkJ+6+q4B/UdVjqtoJfBJ4V1qfJPBPqhpT1RHgz4BPqOpOdXhcVbuAN+C4vL6lqnFV/SPwI+BPMgchIquBc4F/cM97H84P9HT5c+DrqvqQqiZU9VYgBlyQ1uZLqnrAvYc/Bbaq6lZVTarqXcA2HCPk8S1Vfc5tfwewMe3f625V/b6qjqlql6puFxFxx/FXqtqtqgPA/wKumMF9GfMY8+MaBrxJVe/O2Lcc2Je2vc/d59GpqtG07VU4KiCTNcD5ItKbti8AfCdL2+VAj06ck9nnnns6rAGuEZEPpu0LMfE+DmS0f5uIvDFtXxD4Tdr20bTPw4wHNuS6/0VALfCoY38ARwn6i7wHo8owo2MY2TmM8yP8lLu92t3nkZme/QBwIvBklv2/VdVLirjmEaBFROrSDM/qLNcqlgPAp1X103napJ/7APAdVf3zaV7rvCz7jwMjOO7EQ9M4r1FlmHvNMLLzfeATIrLInSz/R+A/8rS/GfhXEVnvrpfZICJtwM+Bk0XkXSISdF/nishpmSdQ1X047qxPikhIRF4OvDGz3RT4BnCdiJzvjqlORF4vIg052v8H8EYReZ0bGBFxgyZWFnGt7wIXi8jb3YCENhHZqKpJdxz/JiKLAURkhYi8bgb3ZcxjzOgYRnY+hWMAduAECPzR3ZeLz+PMcfwK6Ae+CdS4cxivxZnDOIzjnvICELLxTpyJ/27gn4DbpnsDqroNZz7ly0APsAsn2CBX+wPAZuBjQCeOevk7ividUGcN0eXA37hj344TTAHwEffaD4pIP3A3TsCFsQARVSviZhiGYVQGUzqGYRhGxTCjYxiGYVQMMzqGYRhGxTCjYxiGYVQMW6dTgPb2dl27du1sD8MwDGNe8eijjx5X1UWZ+83oFGDt2rVs27ZttodhGIYxrxCRfdn2m3vNMAzDqBhmdAzDMIyKYUbHMAzDqBhmdAzDMIyKYUbHMAzDqBhmdAzDMIyKYUbHMAzDqBhmdIyKsK9riJ89frhwQ8MwqhpbHGqUnUO9I7zj6w9ybCDK689chs8nhTsZhlGVzKrSEZFLRWSniOwSkRuyHBcR+ZJ7fIeInJ127BYROSYiT2b0aRWRu0Tkefe9Je3YR91z7bTKhZWhZ2iUq7/5EEf7oyQVRhPJ2R6SYRizyKwZHRHxA18BLgNOB64UkdMzml0GrHdf1wJfSzv2beDSLKe+AbhHVdcD97jbuOe+AjjD7fdVdwxGmRgejfOebz/CgZ4RLj1jKQCxuBkdw1jIzKbSOQ/Ypap7VHUUuB2nVG46m4Hb1OFBoFlElgGo6n04ZXEz2Qzc6n6+FXhT2v7bVTWmqi/glM89r5Q3ZEzky7/exeMHe/nSFRt52fp2AEbN6BjGgmY2jc4KnBrsHgfdfVNtk8kSVT0C4L4vnuq5RORaEdkmIts6OzsLXM7IxZG+KKtba7n0RcsI+50/tVg8McujMgxjNplNo5NtNlmn0aaU13N2qt6kqptUddOiRZMycxtFEosnCAecP7GQ+25KxzAWNrNpdA4Cq9K2VwKZMbXFtMmkw3PBue/HZnAuYwZEx5KEA860mWd8bE7HMBY2s2l0HgHWi8g6EQnhTPJvyWizBbjajWK7AOjzXGd52AJc436+Bvhp2v4rRCQsIutwghMeLsWNGNkxpWMYRiaztk5HVeMicj1wJ+AHblHVp0TkOvf4jcBW4HKcSf9h4D1efxH5PvAqoF1EDgL/pKrfBD4D3CEi7wP2A29zz/eUiNwBPA3EgQ+oqk0wlJHYWJJw0DE2nuIxpWMYC5tZXRyqqltxDEv6vhvTPivwgRx9r8yxvwu4KMexTwOfnu54jakRiydprAkCpnQMw3CwNDhG2Uh3r43P6Zi4NIyFjBkdo2zE4kmb0zEMYwJmdIyyEbPoNcMwMjCjY5SNWDyRCiQwpWMYBpjRMcpIunttPHrN5nQMYyFjRscoG47RcYxNyNxrhmFgRscoE/FEkkRSs0SvmdExjIWMGR2jLHjGJTWn47c5HcMwzOgYZcIzOpGg417z+YSQ32dKxzAWOGZ0jLLgBQx4bjVw5nVM6RjGwsaMjlEWYmOuey0wXpw1HPBZ9JphLHDM6BhlITWnY0rHMIw0zOgYZSHlXguO/4k5SseMjmEsZApmmRYRH3AWsBwYAZ5S1Y5yD8yY34wrnXH3mikdwzByGh0RORH4CHAx8DzQCUSAk0VkGPg6cKuq2q+IMYno2ORAgnDAb3M6hrHAyad0PgV8DXi/W9cmhYgsBt4JvAu4tXzDM+Yr2QIJQgEfowl7RjGMhUzOOR1VvVJV78s0OO6xY6r6BVWdkcERkUtFZKeI7BKRG7IcFxH5knt8h4icXaiviPyniGx3X3tFZLu7f62IjKQduzHzekbpyFwcCu6czpgZHcNYyBQzp+MHXg+sTW+vqp+fyYXd834FuAQ4CDwiIltU9em0ZpcB693X+TjK6/x8fVX1HWnX+BzQl3a+3aq6cSbjNooj1zqdwVh8toZkGMYcoJhy1T8DosATQCkfU88DdqnqHgARuR3YDKQbnc3Aba7aelBEmkVkGY4BzNtXRAR4O/CaEo7ZKJJsgQSmdAzDKMborFTVDWW49grgQNr2QRw1U6jNiiL7vgLoUNXn0/atE5HHgH7gE6r6u+kP38hHLEsgQSjgtzkdw1jgFLNO55ci8toyXFuy7MucP8rVppi+VwLfT9s+AqxW1RcDfw18T0Qasw5M5FoR2SYi2zo7O7MO3shP7jkdi14zjIVMMUbnQeDH7iR8v4gMiEh/Ca59EFiVtr0SOFxkm7x9RSQAvAX4T2+fqsZUtcv9/CiwGzg528BU9SZV3aSqmxYtWjTF2zJg3Oh42aXBotcMwyjO6HwOeAlQq6qNqtqgqlkVwhR5BFgvIutEJARcAWzJaLMFuNqNYrsA6FPVI0X0vRh4VlUPejtEZJEbgICInIATnLCnBPdhZCEWTxDwCQG/ZSQwDGOcYuZ0ngeezBY6PRNUNS4i1wN3An7gFlV9SkSuc4/fCGwFLgd2AcPAe/L1TTv9FUx0rQFcCPyLiMSBBHCdqnaX8p6McWJjyQnzOeAoHTM6hrGwKcboHAHuFZFfAjFv50xDpt1zbMUxLOn7bkz7rMAHiu2bduzdWfb9CPjRDIZrTIFYPEk46J+wLxzwMxpPoqo4wYWGYSw0ijE6L7ivkPsyjILE4gkiGUrHUz6jieSEUGrDMBYOBY2Oqn4SQEQanE0dLPuojHlPdqUzXrLajI5hLEwKBhKIyIvctS1PAk+JyKMickb5h2bMZ3LN6QA2r2MYC5hiotduAv5aVdeo6hrgb4BvlHdYxnwnFk9MMjrpSscwjIVJMUanTlV/422o6r1AXdlGZFQFsSwuNFM6hmEUE0iwR0T+AfiOu/2nOIEFhpGTWDxJbWhy9BqY0jGMhUwxSue9wCLgv9xXO+56GcPIRTb3mpedwAq5GcbCJa/ScVfw/0BVL67QeIwqITo22b3m5WEzpWMYC5e8SkdVE8CwiDRVaDxGlZBf6ZjRMYyFSjFzOlHgCRG5Cxjydqrqh8o2KmPeExtLTsgwDaTW7ZjSMYyFSzFG5xfuK52S5mEzqo+s0Ws2p2MYC55ijE6zqn4xfYeIfLhM4zGqhKzrdILmXjOMhU4x0WvXZNn37hKPw6giVNVVOjanYxjGRHIqHRG5EngnTonn9Fo1jUBXuQdmzF/GEooqk3OvWfSaYSx48rnXHsApa9COU8jNYwDYUc5BGfMbb85mknvN73ePm9ExjIVKTqOjqvuAfSJyMTCiqkkRORk4FXiiUgM05h+eUcmV8NOUjmEsXIqZ07kPiIjICuAenGwE3y7FxUXkUhHZKSK7ROSGLMdFRL7kHt8hImcX6isi/ywih0Rku/u6PO3YR932O0XkdaW4B2My40YnV+61ykavxRNJjg/GCjc0DKPsFGN0RFWHgbcA/66qbwZOn+mF3WwHXwEuc893pYhknvcyYL37uhb4WpF9/01VN7qvrW6f03HKWJ8BXAp81T2PUWJiY657LWOdjt8nBHxScaXzX388xCv/z28YGbVQbcOYbYoyOiLyEuAqxtfrFBNqXYjzgF2qukdVR4Hbgc0ZbTYDt6nDg0CziCwrsm8mm4HbVTWmqi8Au9zzGCUml9Jx9vkqPqdzoGeYodEEPcOjFb2uYRiTKcbo/CXwUeDHqvqUiJwA/CZ/l6JYARxI2z7o7iumTaG+17vuuFtEpGUK1zNKQMroBCf/eYUCvoornf6RMec9OlbR6xqGMZmCRkdVf6uq/0NVP+tu7ylRChzJdrki2+Tr+zXgRGAjTvSdF3lXzPWchiLXisg2EdnW2dmZrYmRh5R7LTD5zysc8Fd8TmcgGgegfyRe0esahjGZgm4yEfkZk3+c+4BtwNdVNTrNax8EVqVtrwQOF9kmlKuvqnakjf0bwM+ncD3cc9yEUzGVTZs2WcqfKZLPvTYrSsdVOJ7iMQxj9ijGvbYHGMQpUf0NoB/oAE5mZmWrHwHWi8g6EQnhTPJvyWizBbjajWK7AOhT1SP5+rpzPh5vBp5MO9cVIhIWkXU4wQkPz2D8Rg5yhUx7+yo9p+MpHHOvGcbsU0xAwItV9cK07Z+JyH2qeqGIPDXdC6tqXESuB+4E/MAt7pzRde7xG4GtwOU4k/7DuMXjcvV1T/1/RGQjjjrbC7zf7fOUiNwBPA3EgQ+4pRuMEuO5zyJzZU7HNTZ9pnQMY9YpxugsEpHVqrofQERW42QpAJhROJAbzrw1Y9+NaZ8V+ECxfd3978pzvU8Dn57ueI3iiI7Nrei1VCCBzekYxqxTjNH5G+D3IrIbZzJ+HfAXIlIH3FrOwRnzk1xpcGC2lI651wxjrlDQ6KjqVhFZj5P+RoBn04IHvlDGsRnzlFhepeOnt4LrZeKJJIMxL3rNjI5hzDbFLvI8B1jrtt8gIqjqbWUblTGvKbROp5LuNc/ggCkdw5gLFBMy/R2cdS/bAW/iXQEzOkZWPPeaVz8nnXCF3Wvp8zg2pzN/SCSVsUSSSNAyVVUbxSidTcDp7qS+YRQkFk8S8vvw+Savx6200vHUjU9M6cwnPnT7Y9z9dAevWN/OJacv4XVnLKW5NjTbwzJKQDHrdJ4ElpZ7IEb1EBubXDXUw8lIUEml4xiaZU01FjI9j9h9bJCW2hDPHBngIz96gj/95kOzPSSjRBRjdNqBp0XkThHZ4r3KPTBj/hKLJ7LO54DnXqvc8ihP3axoqckbSDAUi3P5F3/HjoO9FRqZkY+e4VEuPLmd33/k1bx900oO9ozM9pCMElGMe+2fyz0Io7qIxZNZI9eg8ut0vHmclS01PLK3m2RSs7r9jvSN8PSRfrbt7WHDyuaKjc+YjKrSMzxGS20IEWGpq1ITScWf5f/OmF8Um/Az9cJZzf/28g/NmK84Rif7n1Yo4GM0kaRSU4Se0lnZUosqDI5mDybwFrRasbfZZ2QswWg8mZrDaa4Jomoh7+Xg4z9+gr+8/bGKXrOokGk3rcw7cYzNC8CPyjgmY54TG0ukqoRmEg74UIWxhBIKlP+ptT8aRwRWNEec7ZExGiPBSe2ibmbsrkGruTPbdA85/wetdc7/U4v73jM8SkudBROUil8/28F3H9pPe324otfNaXRE5GScRJpXAl3Af+JUEX11hcZmzFNi8dyhrp4xGk0kcxqmUtI/MkZ9OEBTTcjdjkPL5HYjrtExpTP79A47iialdNz3nmFTOqViMBbnEz92ciEfH4w587A5XOKlJt+3/lngIuCNqvpyVf13xtfpGEZOnD/g3NFrMF5zp9z0Rx1l01gTSG1nI+VeGzKlM9t4FV5bXGPjvVcyk0W18//u3MmR/ihXnudUeznWX7mHrXxG563AUeA3IvINEbmI7IXQDGMCsXiScBFKpxL0j8RprAmmXGq5wqY999rxAVM6s42naFpqgxPeTemUhj/u7+HWP+zl6gvWcPmZTiWYw72Viw7MaXRU9ceq+g6cnGv3An8FLBGRr4nIays0PmMekn+dji/VphI4SidAU43zw5VrMtpzr3UNxSoW5GBkp8dVm978TbMpnZLyuV/tZElDhL+79FSWNdUAcKRvurU4p04x0WtDqvpdVX0DTrXN7cAN5R6YMX/J516rvNIZm6B0vIzTmXjuvuhYkqFR8yLPJp57rdl9UGiMBPD7JLXfmBl7jw/z0pPaqA8HWNbkBNjMCaMjIttE5IsicqmIRABUtVtVv66qr6nYCI15R/51Ot6cTmWMzkA0TmMkSH3EndMpoHQAuiyYYFbpHR6jIRIg4ObuExGaa4LmXisByaRybCDKkkbH2NSFAzRGAhzpmwPuNeAC4MfAq4DfishWEfmwG9VmGDlx5nQKKZ0KBRKMjNFY4zwpN0QCBQMJYHoRbD1Do+w6NjjtcRrj9AyPpoIHPJprg+ZeKwE9w6OMJZQlDeNh0suba+aG0lHVuKreq6o3qOr5wPuAAeBTIvKYiHx1phd3VdROEdklIpNcduLwJff4DhE5u1BfEfm/IvKs2/7HItLs7l8rIiMist193Zh5PaM0RMfyRa9Vbk4nkVQGYvGUa60xEsyZaTqapnSOT3GtTu/wKG+98QGuueXh6Q+2SijFfJiTjWDiWqqW2hA9Q6Z0ZkqHG6XmKR2ApU2RuaF0RORPPLcagKoeUdVbVPXtOPV1vjuTC4uIH/gKcBlwOnCliJye0ewyYL37uhb4WhF97wJepKobgOeAj6adb7eqbnRf181k/EZu8rnXPKUTq8CczqA7f9PozQ3UBHMqnZEJRqd4pRMdS/Dnt21jT+cQnQMLOwjhN88eY8M//4q+GbrBeoYmLwJtrg3ZnE4JODbgKJrFaUZnWVMNR3rngNIBrgL2i8htInKZ+0MPgKomVfX+GV77PGCXqu5R1VHgdmBzRpvNwG3q8CDQLCLL8vVV1V+pqvc4+yBO8INRIVSV0TxpcLz9laip4xmYRnc+pzESyBMynUxFuBWblSCZVP72B4/zyN4ezlnTwmgiOcFNt9C4d+cxBmJx9nYNzeg82dxrLbXB1KJRY/ocSymdNPdaU4SuodEJar+c5HOvvRk4CbgH+BBwwA2XvrBE114BHEjbPujuK6ZNMX0B3gv8Mm17nesa/K2IvCLXwETkWjeQYltnZ2fhOzFS5KsaCmnutQoYHc/ANETSlE4OoxMbS9DghlYXq3S+9tvd/HzHEW647FTecvaKCddciOw41AdAR//Mnpp7h8doznSv1ZnSKQXe/82itDmdpW4E20z/34olb8i0qvar6q2qehlwJk649L+LyIF8/Yok20LTTN9ErjYF+4rIx3GSk3puwCPAalV9MfDXwPdEpDHbwFT1JlXdpKqbFi1alOcWjExSRqdA9FpFlU6Np3SCDOQImR4ZS1AT9NNWHypa6dz51FHOXdvC+y88IaWSFqrRGUskefpwPzCzH6/ReJLBWDxrIEEsnmTEwtlnRMdAlJba4ITv5/JmZ63O4Qq52IpKfiUiLcBbgHcArZQm4edBYFXa9krgcJFt8vYVkWuANwBXeRVPVTWmql3u50eB3YBF4pUYr1R1oXU6sQrU1PEMTCqQoCaQU+lExxJEgn7a68N0FqF0VJUXjg9x6tJGRGTBG52dRwdSDxwdM0ip0jsycWGoR0sq/5qpnZnQ0R+bEEQA40rnaH9lggnyBRI0iMi7RGQr8AxwLvApHLXwlyW49iPAehFZJyIhnOSimcXhtgBXu1FsFwB9qnokX18RuRT4CPA/VHU47X4WefNSInICTnDCnhLch5GGF5U2J+Z0XAPgGYSmmiADsTiJ5OTJfk/ptNeHilqn0z00ykA0ztr2ugnXWKhGZ8dBx7UW8vs4OgOl40WoTY5eG880bUyfY/3RCUEEAMubKqt08pU2eAG4Eydi7L9VtaTfJlWNi8j17jX8wC2q+pSIXOcevxHYClwO7AKGgffk6+ue+stAGLhLRAAedCPVLgT+RUTiOIlLr1PV7lLek5E+p1Mgeq0i7rUMpeO+D0bjNGX8qEXHkjREArTXh7l/sKvgub3J8nXttYAZnScO9dIYCbBuUf2M3GuZyT49xlPhLMx/31JxbCDGyUsaJuyrCflprg1ytEJrdfIZndWeUhCRGhE5QVV3lvLiqroVx7Ck77sx7bMCHyi2r7v/pBztf4TVASo7Bd1r/sorHS8bgRc63R8dy2J0EixqCNNWF6ZvZIzReP7SCy8cd0T02rb5qXTGEkmC/tKVlnj8QB8bVjZTF/bzwvHpR695C0AnBRKYe23GONkIYixunFw/Z2lj5dbq5Ite8wzOG3ECCP7b3d4oIpluMMMA0gMJsv9pBfw+/D6pyJxOf3SMhnAgVeLYC53OZhiinnutwflx6y5Q4mDv8SH8PmFVq6N0GgpksZ5L9AyNsvGTv+LPbt2WWrcxE6JjCXZ2DLBhZRNLGyMzmtMZzzA9OWQ6/bgxdbqGRkkkddKcDjjBBHMpkOCfcdbF9AKo6nZgbbkGZMxvvDmdXEXcwDFIlVE68ZS6gTSlk9XoJIkEfbTVOU+BhcKmX+gaYmVLTUoteGl2+ubBk/jeriGGRhPc/UwHr/23+/jp9kMzWtT69JF+Ekllw8pmFjdG6BsZm/aaj4LuNat3NG08t+fihslGZ1lTZEZzcVOhGKMTV9W+so/EqAoKudfAmdepzJyOkzjSYzzTdBajE3eUziJX6RQyOnuPD6Vcax5NNcF5oXS8H58vXfli1rbV8eHbt/ODbQenfb4dB3oBOGtVU+operrzOj1Do0SCPmpCEx9aQgEfdSG/KZ0Z4KnaJVnca8uaInRXaIFoMUbnSRF5J+AXkfUi8u/AA2UelzFPKbROxzlWKaUzlqF0vEzTk9fqjIw6IdOe0sm3VkdV2Xt8iHXtpTc6yaRy6Rfu4zsP7pvRefLhub9eemIbP7zuJaxureWuZzqmfb4dB/tY1BBmaWOEpa7Rme6ktJN3LZT1WHNtKBVSbUydbHnXPCpZV6cYo/NB4AwgBnwP6AP+soxjMuYxhTISQCWVzniyTxif7M9UOqqaqnba3lDYvdY5EGNoNMHattoJ+0thdJ4/NsizRwe488mj0+r/5KE+/u4HjxPPk9uuoz9KwCe01oYI+H2ct66VbXu7p+1i23Gojw0rmhCR1FN0xzQrsPYOj6ZcaZm01FkqnJngpcBJz0bgMV5Xp/zBBMUUcRtW1Y+r6rnu6xOqWrnscMa8wiuGls+9Fg74K6h0xt1rdaEAPpk8p+MZwJqgn7qQn3DAR1eeuQMvOmttGZTOtn1OFP9j+3uyricqxNYnjvCDRw/ylJsdIBsd/TEWN4TxuQEW565toWd4jN2dU486G4zF2d05yIaVzQAs8VKqzEDptNYFsx5rsaSfM6JjIEpbXShr1OIyNytBJRJ/5lscepOInJnjWJ2IvFdErirf0Iz5SDHutZDfV7HotXSl4/MJDZHgpOqhXmqVSNCHiNBeH+Z4nid1b43OCe31E/Y7Rid7mp1i2ba3B4Ch0QTPHs1tOHKxv9sJ5X5kb+4laB390ZRxANi0ttW99tSXrT1xsA9V2LCqCYCGcICaoH9Gczq5lE5zbciUzgzItjDUY1kqK8Hsute+CvyDiDwjIj8Qka+KyC0i8jucOZ0G4IdlH6ExryjGvRYOlt+9lkwqg7F4Kkzao7FmcqbpaNwzOo6hbK8PcTyv0hkm6BeWN0/8Aje5CUVnEgm2bV83L1rhpAT8476eKff3jM7DLxQwOmkRTCe019FaF+KRvVO/3pNuks8NKxyjIyIsnUEklJNhOpfSCZrSmQFOCpzJrjVw/vZbaoMc7nXca9GxBL9//nhZxpFvnc52t3bOuTi1a36Hk2rmz1T1LFX9oqpaXV9jAkVFr/nLb3QGYnFUmRBIAF4ht4lGx1M6NSmjU0DpHB9iVWttqpxy6tw1wRmVN+joj3Kge4Q3bVzB4oYwj07D6OzrcozOtn09OY1fR390wo+PiLBpTUvKtTcVDvYM0xAJ0FY/fr7FDeHU/MFUSCaVvpH8gQR9I2PTcjsakx82MlnW5FQQfWRvN5d/6Xe8+1sPp4xQKSlmTmfQrSD6fVX9SamzEhjVhfeDG8qz2j1Ugei1gVQtnSxGJyOQIJpaW+SMua0+RNdQfvfauoxwaZh5VgLPtXbOmhbOWdPCtikanb7hMfpGxjhxUR3dQ6Ps7pxcPntkNEF/ND7JzXLu2lb2dQ1zbIoKpaM/lopY85iu0umPjpFUcgcS1AZRzb7OyshPIqkcH8ytdMBxsT2w+zhvu/EPxMaSfPPd56YyUJeS0uXBMBYsN/9uD1sed5J8x+JOqWo3711WwgF/2ZWOFxadHkjgbWeGTE92r4XpGhwlmeWJOplU9nYNTQoigPHULdM2Ovu6iQR9nLG8iXPWtHCwZ2RKcyOea+1PznESsD/8wmSj5Z0v01BsWtvijmFqhu5of3RSCO6Sxggd/dEpuxm9NTj5AgmcduZimypdgzGSCotyzOkAnLS4nlg8yXtetpZf/dWFvPLk8pR1MaNjzIgnD/XxqV88w9//8HEO9Y4QG8tdNdTDWadT3kCC/qkondGJRqetPkw8qVkXkXYMRImOJbManZkqnUf39XDWymZCAR/nrGlJ7SuWfd1OgMMrT15Ee30oa2CAZ3QyDcUZy5uIBH15AxCycSyH0YnFk1P+d/BSD+UOJJj/qXCSSeW+5zr5q//czi92HKnYdVNrdLKES3t86KL1/O7vX80/vfEM6sL50nLOjKLPLCJ1qjqzOrRG1fF/79xJU02QWDzB/9r6DI2RYM4M0x7hCqzT8VwwmXM6TVmqh2YLJABnrU7mD6AXLl1q99pQLM5Th/v5n688EXCMQDjg49F9PVx+5rKizuHN56xuq2XTmlYezmZ0BiaXKwbH5blxVXPKxVcMXgLJzHOl1ur0T/73y0dvjhQ4Hi2pTNPzT+ns7xrmp9sPccejBzjQ7cyTDI/Gef2G4v5vZ0quh4106sKBshobj4JKR0ReKiJP49TUQUTOEpGvln1kxpznD7u7+O1znXzg1Sdy3StP5Bc7jvDQnq6CSqcSczqZZQ08GmuCDI0mJiye9OZ00gMJAI5nyUqw18su3V476dhMjM7jB3pJJJVzXDdXKODjrJXNU3J37e8apr0+RH04wLnrWjnYMzJpsZ83Z5MtdPbcta08dbiPwVhxYd/Hh2LEk5oqAuaxdJqpcMaTfRZyr80PpbOva4ibf7eHzV+5nwv/72/43F3PsaK5hi9esZEzVzQxXMEqqB0DhY1OpSjGrP0b8DrcImmq+riIXFjWURlzHlXls//9LMuaIlz9krWowh2PHGDP8SFOXDRZBaRTWaWTMafjhlD3R+O0utUp09fpQLrRmRxMsLdriFDAlyp8lc5MjM4je3sQgbNXt6T2nb2mhZt/tydV1bQQ+7uHWe1mvT7PXXvz8AvdbN64ItXmaF+USNA3KZQcnPU6SYXt+3t5+fr2gtc7liOtirc91WCC8bIGOdxr7lzPXFU6w6NxHtnbw/27jnPPMx2pxbanL2vkhstO5Q0blrGyxfn/+eGjBxkq0riXgo7+GCLjKn42KUpLqeqBjIlhK1S+wPnV0x1sP9DLZ996ZuoH8aOXn8YHv/9Y3oWhUHqlk0gqw6NxRkYTDI8mGBlL8PwxJ3KrPpwZSDBuGDyjk+lea3O/mNnyr71wfIg1rbWp1fzpzKS8wbZ93ZyypCFluAA2rWnhxt8qOw72cd661oLn2N89nGp32rIG6kJ+Htk70eh0DDjRZtkCPc5e3YxP4GM/foJlTRFCAR9XnLs6pwvIy6+WaXS8ei1TjYTrGR7F75OsBhGchacBn8yJQIJkUtlzfIgnD/XxxKE+Hj/Qy/YDvcSTStAvXHBCG396wRpec+pi1mRxxdaG/HROM1XQdDjWH6W9PjwpzH82KMboHBCRlwLqlob+EK6rbaa4paW/iFP982ZV/UzGcXGPX45TOfTdqvrHfH1FpBX4T5zyC3uBt6tqj3vso8D7cIzmh1T1zlLcx3xGVemPxukbHqN3ZJReN+zWe/VHx+gfibvvYwzG4gxG4xzuHeHERXW89eyVqXO9YcMy7th2oOBTuRO9NrXnlm/ct4ftB3pT4xiIxhlwxzKSIzNue31o8loa1zAMpAUJjGQEErTUhvBJdqVzoHuYNW2TXWswXt4gV0jvcx0D/O0PHufmazZNSC+fSCqP7e9l88blE9qfvcaLKOueZHSSSWU0kUyNORZPcLhvJFXfJ+D3cfaaFh7JiGDryLMqvSES5PrXrOex/T2MxpPs6Rzi73/4OC85sS1loCecayB7JFw44Cw0nKrS6R4ao6U2mDPyUURorg1W3L0WHUvwfMcgzxzp5+kj/Tx5qI9njvQzNDq+Ju2M5Y38+YUn8NIT2zhnTQu1ofw/rbWhAEOjlVQ6URbnCSKoJMUYnetwftxXAAeBX5GjmudUEBE/zqLTS9zzPiIiW1T16bRmlwHr3df5OKWzzy/Q9wbgHlX9jIjc4G5/REROB67ASV66HLhbRE5W1apSbcmk0jM8yvHBUboGYxwfct67h5x9PUOjdA+P0j3kfO4tsNguFPDRGAnSWBOgIRKkMRJgaWOEs1e38O6XrZ3woy4i3PLuc/HlCZf2zplUiCeSRT15qSr/+5fP0FIbYk1bLc21IVa11tIQCVAfDlAbct/DfmqCfmpDfsJBPydkiTAbr6kz/oX3XH2ee83vE1rrQlnndI4Pxnhxmgssk3z5177ym13sONjHIy/0TFAPe7uGGIzF2biqeUL71roQL1rRyDfu28Mlpy1hvVtmuHtolKtveQiAn13/ckSEQz0jqMKa1nGDeO7aVj5/13MTkmge649y5sqJ10nnry85OfX5+Y4BXvuF+/jG7/bwkUtPndS2oy+KL4fLZsk0irnlS/bp4aTCKa3SUVW6hkY52hd1Xv1RDnQPs7tzkD2dQ+zrHk59R2pDfk5f1sifnLOSM1Y0sWFlEyctqp+ygqgN+RmOle6nJxZP0DfsPCQOjyaIxZNExxKMxpOMJZKu27u+8IkqQEGjo6rHgXLkWDsP2KWqewBE5HZgM5BudDYDt7llqx8UkWYRWYajYnL13Qy8yu1/K3Av8BF3/+1uFoUXRGSXO4Y/lOHeSs5oPMmxgSgd/TGO9Ufp6I9ybCBG50CMzkH3fSCWqg6YiU+cJ/jWOue1fnE9zbUhWuuCtNSGaKpx3ptrgzTVOK/GmmBRcwnpFFMC2Qs0iMWT+EQYiMUnuJUm3XsiSVLhvS9fxwdenbUaedGkyhukKZ3oWAKfTFzQ2lYXpitD6SSSSvfQKIvy+MVzGZ0jfSOpENmdHQO8nnGjs/PoAACnLm2c1O/LV57N277+B666+SF+eN1LqQn5+dObH2Jnh9Nnd+cQJy2uZ5+7RiddhXnq6NF9PVx02hJUlaP9US4u8ol3/ZIG3rBhObc+sJc/e/m6CVkHwJknyOWy8dbqTIV8KXA8WmqD9AxNXemMjCbY3z3M3q4h9nUN8cLxYfZ1DXGod4QjfdFJ7t6Q38fa9lpOWdrA6zcs47RljZy2rJHVrbWpSrQzoS5cWOkcH4yx9YkjqUW/ntdhIOYo/UFX7Q9Ex4rKgvHa05fMeNyloKDREZFbgQ+raq+73QJ8TlXfO8NrrwAOpG0fxFEzhdqsKNB3iaoeAVDVIyKyOO1cD2Y51yRE5FrgWoDVq1cXeTvTZ3g0zuFe5ynrSN+I894fTT15dfRHs2Y99vuERfVhFjWEWdIY4UXLm2hvCLGoPkxbfZj2+jCLGkK01oVprglmnYeYDUKu0fng9x/jsf099I6M8du/fTWrc7itCpXAngqpQm4jE91rkaB/glvHyUow8d+8Z3iUpDLpxzedpppg1ifxWx/YR1KVtroQz7lGxmPn0QFEYP2SyU+ia9vr+O6fnc/bv/4H3nnzg4QCPo70Rvl/bzuLv/3B49zzTAcnLa5nf1q4tMfGVc0E/cLDL3Rz0WlL6I/GiY4lpxTB9OGLTuLnOw5z0317+Ojlp004lm1hqMfSxsiUE5b2Do+lAiFy0Vwb4oBrYNNRVToHY+zvGmZ/9/jrQPewk2khY/6kpTbImrY6Nqxs5tIzIixrirC0qYalTU5NoEUN4ZIYl1zUhvxEx5IkkprzOt++fy9f/s2uVPt0j0NrXYg1bXXUhwM0RAKph8TGiKP8I0EfkaCTMT3od16Z9Z9mi2Lcaxs8gwOgqj0i8uISXDvbv3TmI3quNsX0nc71nJ2qNwE3AWzatGlGiZ5Uld7hMQ72jHCgZ5iDPcMc6hnhUG+Uw70jHO4byZo5t7Uu5BTFaopw1qpm93OYxQ0RFjc67211oTljSKaCVzBq59EBTlnawIN7ujnaH81pdLxqhoXW/xSD514bSMs0HY1Pjg5rqw+z42DvhH3eHE9bAaWT+YQ/FIvzvYf2cdmLlpFIKs91TDQ6z3UMsLatLqeqPHlJA7e99zze+Y2HUFVufe95nLeulVt+/wJ3P9PB+195Ivu6hp3Kp2kGMRL0c9bKZh5yk396E/tLmoo3OictbmDzWcu57Q/7+PMLT0hF9oEzT+BFY2WypDFM50CsaBcqOAXEvMwIuWipDbJtb5RbH9jLvpSBGWJ/9/CEp30Rx/Ctbq3llScvYk1bLavb6ljTWsvatjqaCiiqclPnzvkMj8ZTASiZ9I6M0lwb5OGPXZx6UKsGijE6PhFpSZuMby2yXyEOAqvStlcCh4tsE8rTt0NElrkqZxlwbArXKxl/c8fjPHGol0M9I6kJR4+GSIAVzTWsaK7h7DXNLGuqYXlzxHlvqmFxY3jKbq35xOvOWMLj//hammqDPLK328n1lCewIDZWOqVTF/I7NXUmuNeSqTU6Hu31oUnRa952ewGlk1ne4IePHqQ/Gud9r1jHvTs7+dXTRyeEQe88OsDJWVROOhtWNvOzD74cYbyOz8WnL+HLv36e7qHRVLh05iT8eetauem+PQyPxotalZ6ND120ni2PH+brv93Nx19/emp/R39uI7G4MUJS4f7dXQzH4nQNjfLWs1dOKkPt4QWuFFI6y5tr6Bke45+2PEVN0M/q1lrWtNVx4fpFrG6rZVVrLatba1nRXDOnv0O1YWdsw6OJnEZnOJagPhyoKoMDxRmPzwEPiIhXxuBtwKdLcO1HgPUisg44hDPJ/86MNluA6905m/OBPteYdObpuwW4BviM+/7TtP3fE5HP4wQSrAceLsF9ZEVR1rTV8bKT2lnZ4nwJVrXWsLKlNu/8xUJARFJPmhE3vDqfTzqWEdI802s3ZmQlGBlLTCrF0F4fZjAWn2AcPKWTb61DenkDESGRVL51/wu8eHUzZ69u4UhvlKTC7s5BzljeRHQswd6uId5w1vKc5/TIdI9cctoSvnTP8/zm2WPs7x7KGpp77rpWvnrvbh7b35uKJpvqAsETFtVzyelL2PrE0ZTRiY4l6Bkey5m1eIWbKPKaW8a/YnVhP29+8cqs7T2X2aocysnj/ReeyKtOWczy5giL6sN5c/zNZcaVTu6HrcFYfFLIfzVQTCDBbSLyKPBqHBfVWzIizKaFqsZF5HrgTpyw51tU9SkRuc49fiOwFSdcehdOyPR78vV1T/0Z4A4ReR+wH8dI4p77DpxggzjwgXJGrn3+7RvLdeqqwvuxj+YIe3aOlU7pgJd/LS16bSyRMn4ebW6IcNfQaOoH9HgRSie9vEFNyM99z3Wyt2uYv3udE/11ylJH0TzXMcAZy5vYdWyQpMIpbmTaVHjRikaWNIa56+kO9ncPc+H6yQkaz1nTgk+cRaLeE/N0VqWfs6aFO5/qoGswRlt9OLXGJJer7mUntfPZt55JU02QZU01vOmr9/PC8clzMR4He1yjU0Dp1IT8k6L85iO1ruLLt0B0eDSRaldNFGtGnwV6vPYislpV98/04qq6FcewpO+7Me2zkiM8O1tfd38XcFGOPp+mNCrNKBHej32+DAWlVDrgZZqeqHQy3T5esEDXYCxldLoGYwR8Mim1TjrpWQlqQn7u33WccMDHJW7k0Jq2OoJ+YedRZ/GqF7l2ytKpGx0R4aLTlvCDbQcYS2jWObHGSJDTljXy8AvdnLyknsZIIKeLKx9nrmgG4IlDfbzqlMUFVVMo4OMd544H4SxvqmF/V+7UjV6G7EJGp1rwcpwVUjoNORbKzmeKyb32QaADuAv4OfAL990wZkykCKVTyjkdmJxpOjqWTI3Doy0t6afH8cEYbfX5gzcyU+E8dqCXDSubUioj6Pdx4qL6VDDBcx0DTnhujiCKQlx82mLGEk6sS675kPPWtfLYgR4O9oxMO/fWGW41U69SaK4SCblY3VqbCuvOxoHuERrdKKyFQErp5AmbHh6tTvdaMd/iDwOnqOoZqrpBVc9U1Q3lHpixMAgXoXQy09TMFKd6aFr02lhiUiDBoixJP7sGR2mryz8Jn250RuNJnjjUN2kx6clLGlJGZ2fHACcunvriQo+XntieGnu2OR2A89e1Eh1L8sDurmkbncZIkHXtdTzhGh0vBU6xRmdNW20qrDsb+7uHc0YvViMppZNngehQLFEws8F8pJi/9ANAX7kHYixMipnTKbXSaYgEJigdJ5Ag0702Of/a8aHRvOHSMNHoPHOkn9F4khdnzEGcvKSegz0jDMbi7Dw6wKnTcK15RIJ+XrG+HZ+MT95nsslN/jkylkjlRZsOZ65o4omDzk/BsYEY4YBvUkLVXKxuq6VraDRnBusDPcMFgwiqiWKUztBonPrwwpzT2QPcKyK/AFK+BlX9fNlGZSwYUtkJ8gUSlHxOZ2L0WmwsOSmQoDYUoCbon5CV4PhAjBMLLLBLNzqH3MnxbEoHnEwBR/qiqe3p8jevPYWLT1uSM7S2vT7MiYvq2N05NKPU9meuaGLL44c5PhjjaF+UpU3ZE4dmY02r8++2r2uIM5Y3TTiWTCoHe0a4+LS5sWK+EqSi1/IEEgzF4tRWoXutmDva775C7sswSoaIFCx1UI45Ha+mTsDvc9xrocnnbqsPpeZ0nPxcsSkpnR0He1MLfNPxggZ+7pb4nonS8c5XKBDhvHVt7O4cKtodlo0zVzrG4olDfU42ghzh0tnw0vPs7xqeZHSODcQYjScXTBABjK/TyVzD5xGLJxhLaFXO6RQTMv1JsMqhRvmIBP0FQqZLH70GTlaClroQI1lCpsGJYPNS4QyNJoiOJfOGSzvnHjc6j+3v5cWrmye1WdVSSyTo47+fOgrAyTM0OsVw3roWvv/w/klVPqfCGcvdYIKDfQUTh2bizddkCyY44IVLt2R3D1YjIb8Pv08YzuFe8+Z6qjFkupjotZdY5VCjnBRUOiXMvQZp+deiziLOaJaQaYBF9eOZprtSKXDy/2h75Q1eOO6kZslmdHw+Yf3iBgaicRrCAZZPIS3NdLn4tCVc/ZI1vOTEwsXZctEQCXJCex07XKWzdAoGrNHNF7YvSzBBKm/cAlI6IkJtyM9QjkACb+6rEuWjK00x3+Iv4FQO7QKncihglUONklFY6ZTY6KSVNxhLKEnNrqLSM017xqeQew0cF9vvnu8EJs/neHjzOCcvbajIqvqGSJB/2fyiGYckn7myiYf2dE05cSi4YdNZ1uoc6BlGBFYsIKUDzrxOTqXjut3qFmj0Gqp6IGNXVdWgMWaXSNBXMA1OwCclq3o4XrJ6LFUALptBa6sP0T00SjKpqbmdRQWUDniZpscI+IQXZcxfeHi51mYaRFBpzlzRlMrmMFWjs6atNqvSOdA9wpKGSMGKs9VGbdifc05nXOlU379JUSHT6ZVDReRvKVHlUMOAwpVEncWbpfvyjSudsVTUXDb3Wlt9mHhS6RsZS4VOF6t0AE5b1phz9b83jzPTIIJKc+aKcSM6ZaPTWsuRvpFJtWsOuMlKFxp1oUDO6DVPAS1U99p1OKlovMqhGylB5VDD8ChG6ZTKtQZpRidN6WQLJPASe3YNxcbLGhRYHArjRifbfI7HuWtbufzMpVx02uKcbeYiZ6xowvMGTjUSbnVbHUkdz7PmcaBnmJWtC8u1Bk6QQC6l4+Vkq0b3Wt47cstCf0FVy1E51DAAR+nkWyRXcqXjuddG4iljl+387WlZCboGYzRGikszX4zRqQ8H+OpV50x16LNOfTjAuvY69nQOTXmh6Zq0CLYT3NLJsXjCqaW0EJVOOMCxgezVVb0Ag2oMmc77DXKzMC8SEVufY5SNSNCXWouTjVIrnbpQIFVTJ5pyr2Wf0wEnK8HxoVHai6xDkzI6q/IXJJuvbFzVTHv91Gs+rWkdX6vjcahnBNXCJQ2qkdqQP2fCT+8hrLYK53SKMaN7gftFZAuQCj2xjARGqQgH/amsA9mIjiVLUjXUw+cTGiJOVoJ87jXPlXZ8MMbxgRjtRbjWALc09Fjqyb7auOHSU3nPS2OFG2awqCFMTdA/IZjgQM8IsHCyS6fjzOnkcq9Vr9Ip5o4Ouy8fML9mPY15QThQWaUDbnmDaHx84WmWCf+W2iAizhqdrqFR1i/OX93T47x1rZy3rrWk451LLG6MsHgamQ1EhNWttezvHg+b9koaLET3mhO9lt2tPBSL45PSLROYS1hGAmPWiQTzR6/F4pNLD8yURlfppOZ0siidgN9Ha22I40OjHB+M8ZIT2ko6hoXI6rZa9h4f/xk52D1MKOBj8RRLaFcDzjqdRKrKbDqDsTh14cC8rYyaj1nJSCAirSJyl4g8775ndX6LyKUislNEdonIDYX6i8glIvKoiDzhvr8mrc+97rm2u6/5FTZUxUQC/vzRa2OJkq/haIwEGUhXOjmMWlt9iI6+KL3DY0WFSxv5WdtWy/7uYZJJpwbQ/u5hVrbU5K1RVK3Uhv0kkpo1G8fwaLwqI9dg9jIS3ADco6rrgXvc7Qm4kXNfAS4DTgeuFJHTC/Q/DrxRVc8ErgG+k3Haq1R1o/s6NsN7MEpEOOirvNKpCWQEEmQ3am11YZ4/5lT5LJR3zSjM6rY6YvEkx9xy1wutpEE6qUzTWYIJhmKJqlwYCrOXkWAzcKv7+VbgTVnanAfsUtU9qjoK3O72y9lfVR9T1cPu/qeAiIjYL8UcJxLwM5ZQEu7TbybRMiidYgIJwFE63rxDuymdGeNFsP362WP8763P8NzRQVYtwDU6kFZTJ8sC0aHReFUuDIXiAgkmZCQAPsTMMxIsUdUjAKp6JIerawVOATmPg8D5U+j/VuAxVU0Ps/mWiCSAHwGfUtWsv3Iici1wLcDq1auzNTFKiFfILRbPXimxbHM60fzrdGCiujGlM3O8iL6P/fgJAj7hVacs5t0vXTu7g5olUtVDsyqd6nWvFXNX1wFfZDwjwa8oIiOBiNwNLM1y6ONFji2bkzf7o/Dka58BfBZ4bdruq1T1kIg04BiddwG3ZeuvqjcBNwFs2rSpqGsa0ycS8KqHJqnNIibKoXQaawIMxuKpp8xcUULp6qZQhmmjMKtaarn2whNY3BDmTS9esaANeU2e6qFDsQTLm2eWnHWuktPoiMhnVfUjwKunk5FAVS/Oc+4OEVnmqpRlQLb5lYPAqrTtlTih2wA5+4vISuDHwNWqujttPIfc9wER+R6O+y6r0TEqi6cycmWaLpfSAeh0yy7nmshONzQWSDBzfD7hY5efNtvDmBOMVw/NonSq2L2W75t8uYgEgY+W4bpbcCb6cd9/mqXNI8B6EVnnuvWucPvl7C8izcAvgI+q6v3eiUQkICLt7ucg8AbgyVLekDF9xt1rk6N4vHo3pVc6jtE5NhDNGUQA0FbnGJpQwEdDlf4IGLNDbQGlsxCNzn/jRINtEJF+ERlIf5/hdT8DXCIizwOXuNuIyHIR2QqgqnHgeuBOnDmkO1T1qXz93fYnAf+QERodBu4UkR3AduAQ8I0Z3oNRIrxJ/GxKJ5706t2UWuk4X+hjA7GcQQQwrnTa60JVuWbCmD3G53SyGZ04dVVYNRTyz+l8QlX/TkR+qqqb87SbMqraBVyUZf9h4PK07a3A1in0/xTwqRyXnX/ZFRcI+ZRONFXvpjxKp6M/Rn2e0FRvTqfYvGuGUSx1qei1iQ9biaQyMrYwlc4f3PeZqhrDyEs+peMZonLN6XQNxfImrvQmuj03m2GUitocSidVS2cBRq+FROQa4KUi8pbMg6r6X+UblrGQCOcJJCif0nH+9DVHqWqP2pCfSNBnkWtGyakJZlc63na1Kp18d3UdcBXQDLwx45gCZnSMkuCFK2dzr3n7wiXPSDAejppPRYkIH3/96bxoeWNJr28Yfp9QE/SnFih7DI1Wb6lqyGN0VPX3wO9FZJuqfrOCYzIWGPlCpsuldOpDAUQcpVNToGzCuy5YU9JrG4ZHXdg/KSNBNVcNhfzrdF6jqr8Gesy9ZpSTlNLJkvSzXHM6Pp/QEHbKG5SyKqlhTIVaN9N0Op57rRoLuEF+99orgV8z2bUG5l4zSoj3o58t6We5lA44Lrb+aLyg0jGMclEbyq10qrGAG+R3r/2T+/6eyg3HWIh4KiZbeYNyKR3wIthGSlqV1DCmQl04i9LxSlUvQPfaX+fraOWqjVLhqZhsSidWVqXj/PmXw6AZRjHUhvwMTlI61VuqGvK717zS1KcA5zKeguaNwH3lHJSxsAj6BZ/MltIpHEhgGOWiLhTgWH9swr7hBRy95pWp/hVwtqoOuNv/DPygIqMzFgQiQiTozx+9VgbD4IVNWyCBMVvUhvyTcq95yqda3WvFPD6uBkbTtkeBtWUZjbFgCQd8edfpRHKUHpgJDW7+NVM6xmxRG/ZPmtMZHk1QE/Tjr9IS3sWY0u8AD4vIj3Gi1t7MeNVOwygJs6J0Ip7SsTkdY3aoCwUmRa8NxuJV61qDIoyOqn5aRH4JvMLd9R5Vfay8wzIWGjmVzlj5lI7nXrPoNWO2qA0FiMWTxBNJAn7nb3w4Vr21dKA4pYOq/hH4Y5nHYixgciqdeAK/T1JfyFLSaO41Y5bxFM3wWIJG9298MJa9bHu1YH4FY04QDvqJ5lA6uUpJzxQLJDBmm9os1UOHYvG85TbmO2Z0jDlBOOBLrclJJxpPlM0oeLVyPMVjGJXGUzrpEWzDo3FTOiKyRkQudj/XiEhDoT4FztcqIneJyPPue0uOdpeKyE4R2SUiNxTqLyJrRWQkrWrojWl9zhGRJ9xzfUmsDOScIjILSufs1S18693ncu7a1rKc3zAKkU3pDMbiVbswFIowOiLy58APga+7u1YCP5nhdW8A7lHV9cA97nbmdf3AV4DLgNOBK0Xk9CL671bVje7rurT9XwOuBda7r0tneA9GCcmtdJJlUzoiwqtPXYyvSkNTjblPqnroBKWTqOrotWIeIT8AvAy3gqiqPg8snuF1NzMedn0r8KYsbc4DdqnqHlUdBW53+xXbP4WILAMaVfUPqqrAbYX6GJUlEvTniF5LlE3pGMZs41UPHRmdqHQWunst5v7oAyAiAZz1OjNhiaoeAXDfsxmxFcCBtO2D7r5C/deJyGMi8lsR8cK8V7j9s51rEiJyrYhsE5FtnZ2dU7kvY5pEAr4c0WtJC2k2qpZMpaOqDI8mqtq9Vsyd/VZEPgbUiMglwF8APyvUSUTuBpZmOfTxIseWzedRyNgdAVarapeInAP8RETOmOq5VPUm4CaATZs2zdTAGkUQDuZap2NKx6hePKXjzenE4kkSSa3aWjpQnNG5AXgf8ATwfmArcHOhTqp6ca5jItIhIstU9Yjr+jqWpdlBYFXa9krgsPs5a39VjQEx9/OjIrIbONk918oc5zLmAJFArnU6SZrSSksbRjWRqXSqvZYOFOdeqwFuUdW3qeqfALe4+2bCFuAa9/M1wE+ztHkEWC8i60QkBFzBeKbrrP1FZJEbgICInIATMLDHdcENiMgFbtTa1TmuacwS3uJQZ8ptHFM6RjWTil5z53RSVUMX+JzOPUw0MjXA3TO87meAS0TkeeASdxsRWS4iWwFUNQ5cD9wJPAPcoapP5esPXAjsEJHHcSLurlPVbvfY/8RRaLuA3cAvZ3gPRgkJB3wkFeLJDKNTxug1w5htQgEfQb+kFI6neKp5cWgx5jSiqoPehqoOikjtTC6qql3ARVn2HwYuT9veiuPOK7b/j4Af5bjmNuBF0x+1UU48wxIdSxBMS3ljSseodmqC/jSlU91lDaA4pTMkImd7G+4E/Uj5hmQsRMJupufMYAJnnY4ZHaN6qQuPZ5r2auks9ISfHwZ+ICLexPsy4B3lG5KxEIkExpVOOo7SqV5Xg2HUhsaVjvdezYtD8xodd1L+FcCpOGWrBXhWVccqMDZjAeEpncyS1aZ0jGqnLhxIzeWklM5Cda+pagLYrKpjqvqkqj5hBscoB56aicXHlU484axZMKVjVDNNNUGe7xhkKBZn2EKmAbhfRL4sIq8QkbO9V9lHZiwoIlmUjpcA1JSOUc38xatO4kjfCJ/4yZMMue61hb449KXu+7+k7VPgNaUfjrFQyaZ0vASgpnSMauYlJ7bx4YtO5t/ufo4T2usI+qWq/+aLKVf96koMxFjYeGomZkrHWIBc/5qTeOiFLh7Y3VX1GTiKKW3QJCKf9xJgisjnRKSpEoMzFg7p63Q8TOkYCwW/T/jCOzbSXh+iocqLChZzd7cATwJvd7ffBXwLeEu5BmUsPLwFoOnrdLz5HVM6xkJgcWOE//iz8znWH5vtoZSVYozOiar61rTtT4rI9jKNx1igZFU6cVM6xsLi1KWNnJotN38VUcwj5IiIvNzbEJGXYRkJjBLjKZ10o+MpnbApHcOoGopROtcBt6XN4/QwnuHZMEqCp3TS3WumdAyj+shpdERktaruV9XHgbNEpBFAVfsrNjpjwTCudGxOxzCqmXzf5p94H0TkR6rabwbHKBcBv4+ATyau0zGlYxhVRz6jk17i+YRyD8QwnEJuae41UzqGUXXk+zZrjs+GURYiQR9RUzqGUdXkMzpniUi/iAwAG9zP/SIyICIzcrOJSKuI3CUiz7vvLTnaXSoiO0Vkl4jcUKi/iFwlItvTXkkR2egeu9c9l3ds8UzuwSg94YB/YkYCUzqGUXXk/Darql9VG1W1QVUD7mdvu3GG170BuEdV1+OUw74hs4FbVuErwGXA6cCVInJ6vv6q+l1V3aiqG3EWse5V1e1pp73KO66qx2Z4D0aJCZvSMYyqZ7YeITcDt7qfbwXelKXNecAuVd2jqqPA7W6/YvtfCXy/ROM1KkA2peMTCPolTy/DMOYTs2V0lqjqEQD3PZurawVwIG37oLuv2P7vYLLR+ZbrWvsHEcn5SyYi13q55jo7O4u7I2PGRIK+SdFr4YCfPP9VhmHMM8qWWU5E7gayJXT4eLGnyLKvqIAGETkfGFbVJ9N2X6Wqh0SkAfgRjvvttmz9VfUm4CaATZs2WRBFhYgE/JMyEth8jmFUF2UzOqp6ca5jItIhIstU9YiILAOyza8cBFalba8EDrufC/W/ggyVo6qH3PcBEfkejvsuq9ExZodw0MfQUDy17SkdwzCqh9l6jNzCeCqda4CfZmnzCLBeRNaJSAjHkGwp1F9EfMDbcOaAvH0BEWl3PweBN+BkzjbmEKZ0DKP6ma1v9GeAS0TkeeASdxsRWS4iWwFUNQ5cD9wJPAPcoapP5evvciFwUFX3pO0LA3eKyA5gO3AI+EaZ7s2YJuGgb+LiUFM6hlF1zEq1IFXtAi7Ksv8wcHna9lZga7H93WP3Ahdk7BsCzpnRoI2yEwn4JwQSmNIxjOrDvtHGnCFiSscwqh4zOsacIRycrHSslo5hVBf2jTbmDJGAo3RUnSj1WDxpSscwqgwzOsacIZxRyC02lrA5HcOoMuwbbcwZvEJuKaNjSscwqg4zOsacIVWy2l2rEzWlYxhVh32jjTmDKR3DqH7M6BhzBk/pRE3pGEbVYt9oY84wbnSSxBNJ4kk1pWMYVYYZHWPOMO5eS6RcbKZ0DKO6mJU0OIaRDU/p/GF3F995cB8ANSFTOoZRTZjRMeYMntL53F3P0RAO8M7zV/P6M5fN8qgMwyglZnSMOcNpyxp578vWcebKRi49Y5mpHMOoQszoGHOGUMDHP77x9NkehmEYZcRmaQ3DMIyKYUbHMAzDqBizYnREpFVE7hKR5933lhztLhWRnSKyS0RuSNv/NhF5SkSSIrIpo89H3fY7ReR1afvPEZEn3GNfEhEp3x0ahmEY2ZgtpXMDcI+qrgfucbcnICJ+4CvAZcDpwJUi4jn8nwTeAtyX0ed04ArgDOBS4KvueQC+BlwLrHdfl5b4ngzDMIwCzJbR2Qzc6n6+FXhTljbnAbtUdY+qjgK3u/1Q1WdUdWeO896uqjFVfQHYBZwnIsuARlX9gzrFWm7LcU3DMAyjjMyW0VmiqkcA3PfFWdqsAA6kbR909+UjV58V7ueiziUi14rINhHZ1tnZWeCShmEYRrGULWRaRO4GlmY59PFiT5Fln06zz5TOpao3ATcBbNq0qdA1DcMwjCIpm9FR1YtzHRORDhFZpqpHXNfXsSzNDgKr0rZXAocLXDZXn4Pu56mcyzAMwygxs7U4dAtwDfAZ9/2nWdo8AqwXkXXAIZwAgXcWcd7vicjngeU4AQMPq2pCRAZE5ALgIeBq4N+LGeijjz56XET2FdN2DtAOHJ/tQcwCdt8LC7vv+cGabDvFmVevLCLSBtwBrAb2A29T1W4RWQ7crKqXu+0uB74A+IFbVPXT7v434xiNRUAvsF1VX+ce+zjwXiAO/KWq/tLdvwn4NlAD/BL4oM7GzZcREdmmqpsKt6wu7L4XFnbf85tZMTpGeaiWP8qpYve9sLD7nt9YRgLDMAyjYpjRqS5umu0BzBJ23wsLu+95jLnXDMMwjIphSscwDMOoGGZ0DMMwjIphRmceU2y2bretX0QeE5GfV3KM5aCY+xaRVSLyGxF5xs1I/uHZGGspyJVtPe24uJnTd4nIDhE5ezbGWWqKuO+r3PvdISIPiMhZszHOUlPovtPanSsiCRH5k0qOb6aY0ZnfFMzWncaHgWcqMqryU8x9x4G/UdXTgAuAD6RlKZ83FMi27nEZ49nTr8XJqD6vKfK+XwBeqaobgH+lCibai7xvr91ngTsrO8KZY0ZnflNMtm5EZCXweuDmygyr7BS8b1U9oqp/dD8P4BjcQglj5yI5s62nsRm4TR0eBJrd9FLzmYL3raoPqGqPu/kgE1NdzVeK+f8G+CDwI7KnEJvTmNGZ3xSTrRucrA5/DyQrNK5yU+x9AyAia4EX46RAmm8Uk219OhnZ5zpTvaf34WQame8UvG8RWQG8GbixguMqGbOVe80okplm6xaRNwDHVPVREXlVCYdWVkqQpdw7Tz3OE+Ffqmp/KcZWYYrJkD6djOxznaLvSURejWN0Xl7WEVWGYu77C8BH3JyS5R9RiTGjM8cpQbbulwH/w81jFwEaReQ/VPVPyzTkklCC+0ZEgjgG57uq+l9lGmq5KSbb+nQyss91ironEdmA4za+TFW7KjS2clLMfW8CbncNTjtwuYjEVfUnFRnhDDH32vzGy9YNObJ1q+pHVXWlqq7FydT967lucIqg4H2L8438JvCMqn6+gmMrNals6yISwvk/3JLRZgtwtRvFdgHQ57kf5zEF71tEVgP/BbxLVZ+bhTGWg4L3rarrVHWt+53+IfAX88XggBmd+c5ngEtE5HngEncbEVkuIltndWTlpZj7fhnwLuA1IrLdfV0+O8OdPqoaB67HiVJ6BrhDVZ8SketE5Dq32VZgD0559m8AfzErgy0hRd73PwJtwFfd/99tszTcklHkfc9rLA2OYRiGUTFM6RiGYRgVw4yOYRiGUTHM6BiGYRgVw4yOYRiGUTHM6BiGYRgVw4yOYRiGUTHM6BiGYRgV4/8D8IEIwTbUdZcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1)\n",
    "df = d[1:-1,3]-f1\n",
    "plt.plot(x1[:],df[:])\n",
    "plt.ylabel('Force difference (eV/Angstrom)')\n",
    "plt.title(\"Force difference\")\n",
    "#plt.savefig(\"dF.jpg\")\n",
    "#plt.xlim([0.0,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Energy')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEICAYAAACTVrmbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuXElEQVR4nO3dd3ydZfnH8c836U53k86kiy5a2lJ6aJlVsCxlFKmAylKgovxQVGQKAopCUQqKCGUoVKQgGwHBIiK7pCPdMy3dbdKRNM1Ort8f5xRDOEnT9pzz5CTX+/U6r57nftZ1d125x3M/MjOcc865WEoJOgDnnHNNjycX55xzMefJxTnnXMx5cnHOORdznlycc87FnCcX55xzMefJxTnnXMx5cnEuxiStlVQiqajG5/6g43IukVoEHYBzTdQZZjYrXheX1MLMKuN1fecOlrdcnEsQSZdIek/SbyXtlLRG0mk19neS9KikzZI2SvqVpNQa574vaZqkHcCtkrpJekVSoaRPIse/Fzn+j5J+V+v+r0i6OpF1ds2XJxfnEms8sBxIB6YCj0pSZN/jQCUwCBgDnAxcVuvcXKA7cAfwR2AP0BO4OPKhxrW+KSkFQFI68BXgqbjUyrlaPLk4Fx8vStpV43N5pPxTM3vYzKoIJ4BeQA9JPYDTgKvNbI+ZbQOmAefXuOYmM/tDpDusHDgH+IWZFZvZksj1ADCz2UAB4YRC5Dr/MbOtcayzc5/xMRfn4mNS7TEXSZcAW/Zum1lxpNHSHugKtAQ2/68hQwqwvsYlan7PIPzvt679EE42FwD/ivx634FVxbn958nFucZhPVAGpNczUF9zCfM8wl1omcCKSFlWreP/CiySNBo4FHgxZtE6tw/eLeZcI2Bmm4E3gd9J6igpRdIhkr5Ux/FVwPOEB/bbSRoGXFTrmA3AJ8AM4DkzK4lvLZz7H08uzsXHK7Wec3mhAedcBLQClgA7gWcJj8nU5f+AToS72mYQHqwvq3XM48DIyH7nEkb+sjDnmgZJdwE9zeziGmUTCHeP9Tez6sCCc82Ot1ycS1KShkkapbBxwKXACzX2twR+BDziicUlmicX55JXB8LjLnuAZ4DfAS8BSDoU2EW4W+3eYMJzzZl3iznnnIs5b7k455yLOX/OBUhPT7f+/fsHHYZzziWVOXPm5JtZRrR9nlyA/v37k52dHXQYzjmXVCR9Wtc+7xZzzjkXc55cnHPOxZwnF+ecczHnycU551zMeXJxzjkXc55cnHPOxZwnF+ecczHnycU555qpv7y/hndX5sXl2p5cnHOuGVq2pZBfvbqUF+ZujMv1Pbk451wzU1VtXP/cQjq2bcnPTx8el3t4cnHOuWZmxodrmb9+Fzeffihd01rF5R6BJhdJ10gySemR7XGS5kc+OZLOrnHseZIWSFosaWo917xB0ipJyyWdkoh6OOdcsti0q4S731jOhCEZTDq8T9zuE9jClZKygJOAdTWKFwEhM6uU1AvIkfQK4feE3w2MNbM8SY9L+oqZvVXrmsOB84ERQG9glqQhZlaViDo551xjZmbc/OIiqg3umHQYkuJ2ryBbLtOAa4HP3lZmZsVmVhnZbFNj30BghZntndYwCzgnyjXPAmaaWZmZrQFWAePiEbxzziWbfyzYzFvLtvHTk4eQ1bVdXO8VSHKRdCaw0cxyouwbL2kxsBC4IpJsVgHDJPWX1AKYBGRFuXQfYH2N7Q2RsmgxTJGULSk7Ly8+U/Gcc66x2LmnnFtfXsyozE5cckz/uN8vbt1ikmYBPaPsugm4ETg52nlm9jEwIvIO8MclvW5mOyV9H3gaqAY+INya+cJto12yjvtMB6YDhEIhf9ezc65Ju+O1pRSUVDDj0vG0SI1/uyJuycXMJkYrlzQSGEB4PAUgE5graZyZbalx/lJJe4DDgGwzewV4JXKNKUC0cZQNfL5FkwlsikF1nHMuab27Mo9n52zgyhMOYXjvjgm5Z8K7xcxsoZl1N7P+ZtafcEI4wsy2SBoQ6fZCUj9gKLA2st098msX4AfAI1Eu/zJwvqTWkgYAg4HZ8a6Tc841VsXlldz4wkIGpqdx1YmDE3bfxvaa4+OA6yVVEO7++oGZ5Uf23SdpdOT77Wa2Aj4bvwmZ2S1mtljSM8ASoBK40meKOeeas3veXMH6HSU8PeUo2rRMTdh9ZebDDaFQyLKzs4MOwznnYipn/S7OfuB9zh/Xl1+fPTLm15c0x8xC0fb5E/rOOdcElVdWc91zC8jo0JrrTxuW8Ps3tm4x55xzMfDQO6tZtmU3D18UomOblgm/v7dcnHOuiVm1bTd/+PcqvjaqFycN7xFIDJ5cnHOuCamqNq59dgHtWqdy6xkjAovDk4tzzjUhMz5cy9x1u7jl9OFkdGgdWByeXJxzronYsLOYqW8s50tDMjh7TPxWPG4ITy7OOdcEmBk3vrAIAXecHd8VjxvCk4tzzjUBz83dyH9X5HHtqcPI7BLfFY8bwpOLc84luW27S/nlP5ZwZP8uXHhUv6DDATy5OOdc0rvlxcWUVFRx5zmjSEkJtjtsL08uzjmXxF5buJl/Lt7CjycO4ZCM9kGH8xlPLs45l6R27innlpcWMbJPJy4/fkDQ4XyOL//inHNJ6rZXFrOrOHEvANsfjSsa55xzDTJryVZenL+JK08YxKG9EvMCsP3hycU555JMQUkFN724kGE9O3DlCYOCDicq7xZzzrkkc8erS8gvKueRi46kVYvG2UZonFE555yL6p0VeTyTvYEpEwYyMrNT0OHUyZOLc84lid2lFdzw3AIGdW/Pj74yOOhw6uXdYs45lyR+/doythSW8uz3j6FNy9Sgw6lXoC0XSddIMknpke1xkuZHPjmSzq5x7HmSFkhaLGlqHdfrL6mkxjUeTFRdnHMunt5bmc9Ts9dx2fEDOaJvl6DD2afAWi6SsoCTgHU1ihcBITOrlNQLyJH0CtAJuBsYa2Z5kh6X9BUzeyvKpVeb2eHxjt855xKlqKyS655bwMD0NH5y0pCgw2mQIFsu04BrAdtbYGbFZlYZ2WxTY99AYIWZ5UW2ZwHnJCpQ55wL0p2vL2VTQQlTJ49q9N1hewWSXCSdCWw0s5wo+8ZLWgwsBK6IJJtVwLBIt1cLYBKQVcflB0iaJ+kdScfHqQrOOZcQH6zK568frePSYwcQ6t816HAaLG7dYpJmAT2j7LoJuBE4Odp5ZvYxMELSocDjkl43s52Svg88DVQDHxBuzdS2GehrZtsljQVelDTCzAqjxDcFmALQt2/f/a+gc87FWVFZJT97Ntwdds0pQ4MOZ7/ELbmY2cRo5ZJGAgMIj6cAZAJzJY0zsy01zl8qaQ9wGJBtZq8Ar0SuMQWoinLPMqAs8n2OpNXAECA7yrHTgekAoVDIau93zrmg/ea1cHfYs1ccnTTdYXslfEDfzBYC3fduS1pLeBA/X9IAYH1kQL8fMBRYGzmuu5ltk9QF+AFwbu1rS8oAdphZlaSBwGAgN951cs65WHt/VT5PfryOy44bwNh+ydMdtldje87lOOB6SRWEu79+YGb5kX33SRod+X67ma2Az8ZvQmZ2CzABuF1SJeGWzRVmtiOxVXDOuYOzu7SCa5O0O2wvmXmPUCgUsuzsL/ScOedcIG54fgFPf7KeZ79/TKN+pkXSHDMLRdvny78451wj8s6KPJ6avZ4pEw5p1IllXzy5OOdcI1FYWsH1zy1gcPf2XD2xca8dti+NbczFOeeardtfWcK23WU8eMHYpJsdVpu3XJxzrhH415KtPDtnAz/48iGMzuocdDgHzZOLc84FbMeecm54fiGH9urIVScmd3fYXt4t5pxzAbv5pUUUlJQz49JxjfbNkvuradTCOeeS1Ms5m3h1wWaunjiEQ3t1DDqcmPHk4pxzAdlaWMrNLy7i8KzOfG9CtOUSk5cnF+ecC4CZcd1zCyirrOKec0fTIrVp/XfctGrjnHNJYuYn6/nP8jxuOO1QBma0DzqcmPPk4pxzCbZ+RzG/+scSjh3UjQuP6hd0OHHhycU55xKoqtr4yTPzSZGYOnk0KSkKOqS48KnIzjmXQA+/m8sna3dyz7mj6dO5bdDhxI23XJxzLkGWbCrkd28u57TDenL2mD5BhxNXnlyccy4Byiqr+Mkz8+nUthV3nD2SyJt4myzvFnPOuQS4580VLNuym0cvDtE1rVXQ4cSdt1yccy7OPsrdzvR3c/nmuL585dAeQYeTEJ5cnHMujgpLK/jpMzn069qOn3/t0KDDSRjvFnPOuTj6xUuL2VJYyrNXHE1a6+bzX663XJxzLk7+sWATL8zbyP+dMIgxSfzK4gMRaHKRdI0kk5Req7yvpCJJ19QoGytpoaRVkn6vOqZaSLohcsxySafEuw7OORfN5oISbnx+IaOzOvN/Jw4KOpyECyy5SMoCTgLWRdk9DXi9VtmfgCnA4Mjn1CjXHA6cD4yI7H9AUnK/K9Q5l3Sqq42fPpNDZbVx33mH07KJLUrZEEHWeBpwLWA1CyVNAnKBxTXKegEdzexDMzPgCWBSlGueBcw0szIzWwOsAsbFJXrnnKvDI+/l8sHq7fzijOH0T08LOpxABJJcJJ0JbDSznFrlacB1wG21TukDbKixvSFSVlsfYH0DjkPSFEnZkrLz8vL2swbOORfdkk2F3P3Gck4Z0YNzQ1lBhxOYuE1dkDQL6Bll103AjcDJUfbdBkwzs6JaQyrRxlcsSllDj8PMpgPTAUKhUNRjnHNuf5SUV/HDmfPo0q4Vv/n6qCb/FH594pZczGxitHJJI4EBQE7kNz4TmCtpHDAemCxpKtAZqJZUCjwXOW6vTGBTlMtvALIacJxzzsXcr19byqptRcy4dFyzeAq/PgmfdG1mC4Hue7clrQVCZpYPHF+j/FagyMzuj2zvlnQU8DFwEfCHKJd/GfibpHuA3oQH/mfHpybOOfc/by3dyoyPPuWy4wZw/OCMoMMJXDI90fN94C9AW8IzyV6Hz8ZvQmZ2i5ktlvQMsASoBK40s6qA4nXONRPbdpdy7bMLOLRXR3526tCgw2kUFJ581byFQiHLzs4OOgznXBKqrjYu/vNsZq/ZwT+uOo7BPToEHVLCSJpjZqFo+5rf5GvnnIuhR99bw7sr87n59OHNKrHsiycX55w7QIs2FjD1jWWcPLwH3x7fN+hwGpUGJxdJaf60u3POhe0pq+SHT82jW1pr7jqneU87jqbO5CIpRdK3JL0qaRuwDNgsabGkuyUNTlyYzjnXuNz68mLWbN/DPeeNpkszn3YcTX0tl7eBQ4AbgJ5mlmVm3QlPF/4IuFPSBQmI0TnnGpWX5m/k73M2cOWXB3HMIen7PqEZqm8q8kQzq6hdaGY7CD/U+JyklnGLzDnnGqF124u56YVFjO3XhasnegdOXeprubwQ6Rarc9W1aMnHOeeaqoqqaq6aOY8UwX3nH06LZrjacUPV9zvzMHAGsEbS05ImSfKORedcs/XbN5eTs34Xd50ziswu7YIOp1GrM7mY2Utm9k2gH/A8cDGwTtJjkk5KVIDOOdcY/Gf5Nh56J5dvje/LaSN7BR1Oo7fPNp2ZlZjZ02Z2NuGVjMcA/4x7ZM4510hsLSzlJ8/kMKxnB245fXjQ4SSFfSYXST0kXSXpfeBF4E1gbLwDc865xqCq2vjRzHmUlFdx/7eOoE1Lf9yvIeqcLSbpcuCbwFDC3WLXmtn7iQrMOecagz/8eyUf5e7g7smjGNS9fdDhJI36piIfA9wJzDKz6gTF45xzjcYHq/K5762VnD2mD5PHZu77BPeZOpOLmX0HQGEXAAPN7HZJfQk/VOnvSXHONVnbdpfyw5nzGZiexq8mHebLu+ynhkzSfgA4mnAXGcBu4I9xi8g55wJWVW38+On5FJVV8MC3x5LWOplefdU4NOR3bLyZHSFpHoCZ7fTnXZxzTdn9/17F+6u2M/WcUQzt6cvoH4iGtFwqIqshG4CkDMDHYJxzTdL7q/K5960VfH1MH74R8nGWA9WQ5PJ74AWgu6Q7gPeAX8c1KuecC8DWwlJ+NHMegzLa86uzfZzlYOyzW8zMnpQ0B/gKIGCSmS2Ne2TOOZdAlVXVXPW3eewpq2LmlCNo18rHWQ5Gfc+5tDezIgAzW0b4fS51HuOcc8ns7jeXM3vtDu4973AGdfdxloNVX7fYS5J+J2lCzZWRJQ2UdKmkN4BTD+bmkq6RZJLSa5X3lVQk6ZoaZWMlLZS0StLvFaW9Kqm/pBJJ8yOfBw8mPudc8/Dm4i2frRs2aUyfoMNpEup7zuUrkr4KfA84VlJXoAJYDrwKXGxmWw70xpKygJOAdVF2TwNer1X2J2AK4ReVvUY4sdU+BmC1mR1+oHE555qXT7fv4ad/z2Fkn06+blgM1dupaGavEf6PPB6mAdcCL9UslDQJyAX21CjrBXQ0sw8j208Ak4ieXJxzrkFKK6q44q9zSZF44Nu+blgsBfKmG0lnAhvNLKdWeRpwHXBbrVP6ABtqbG+IlEUzQNI8Se9IOr6eGKZIypaUnZeXt/+VcM4lNTPj5hcXsXRzIfeedzhZXf39LLEUt+kQkmYBPaPsugm4kfDy/bXdBkwzs6JaQyrR5gNalLLNQF8z2y5pLPCipBFmVviFk82mA9MBQqFQtGs555qwmZ+s5+9zNnDViYM4YVj3oMNpcuKWXMxsYrRySSOBAUBOJIFkAnMljQPGA5MlTQU6A9WSSoHnIsftlQlsinLPMqAs8n2OpNXAECA7RtVyzjUBOet38YuXFjNhSAZXTxwSdDhN0j6Ti6TfAn82s8WxuKGZLQQ++zFB0logZGb5wPE1ym8Fiszs/sj2bklHAR8DFwF/iBJrBrDDzKokDQQGEx6/cc45AHbsKef7f51DRofW3Hfe4aSm+IOS8dCQMZdlwHRJH0u6QlKneAdVh+8DjwCrgNVEBvMlnSnp9sgxE4AFknKAZ4ErzGxHEME65xqfyqpqfvjUPPL3lPPgBWPpkubLJMaLzBo23CBpKPAdwqsjvw88bGZvxzG2hAmFQpad7T1nzjV1v3ltKQ/9N5epk0dxbigr6HCSnqQ5ZhaKtq9Bs8UiC1cOi3zygRzgJ5JmxixK55yLo1dyNvHQf3O58Kh+nlgSoCFjLvcAZwJvAb+u8ZKwuyQtj2dwzjkXC8u2FHLtswsI9evCzf6gZEI0ZLbYIuDnZlYcZd+4GMfjnHMxtau4nClPzKFDmxY88O0jaNUikMf7mp2GJJf5wLBaz50UAJ+aWUE8gnLOuViorKrmqqfmsaWglKemHEX3jm2CDqnZaEhyeQA4AlhA+GHGwyLfu0m6wszejGN8zjl3wO765zLeXZnPXeeMZGy/LkGH06w0pH24FhhjZiEzGwuMIdxVNhGYGsfYnHPugL04byMPv7uGi47ux3lH9g06nGanIcllWM0HKM1sCeFk4w8nOucapZz1u7juuQWMG9DVB/AD0pBusRWS/gTsnXZ8XqSsNeEl+J1zrtHYVljKlBnZpLdvzZ++fQQtU30APwgN+V2/mPBT8VcDPya8nMolhBPLCfEKzDnn9ldpRRWXz5jD7tJKHrk4RLf2rYMOqdmqt+USeXjylcgilL+Lcoi/4tg51yiYGTc+v5Cc9bt48IKxHNqrY9AhNWv1tlzMrAooDnA9Meeca5AH38nl+Xkb+clJQzj1sGhv+3CJ1JAxl1JgoaR/UePtkGb2w7hF5Zxz++GNxVuY+sYyzhjdm6tOHBR0OI6GJZdXIx/nnGt0Fm0s4OqZ8xmV2Zm7J4+i1gPfLiD7TC5m9riktoTf8OhriTnnGo1thaVc/kQ2ndu15OELx9KmZWrQIbmIfc4Wk3QG4SVg/hnZPlzSy3GOyznn6lVSXsVlT2RTUFLBwxeFfGmXRqYhU5FvJbxA5S4AM5tP+DXFzjkXiOpq48dPz2fhxgJ+f/4YDuvjc44am4Ykl8ooC1Q27A1jzjkXB1PfWM4/F2/h518bzsThPYIOx0XRoCX3JX0LSJU0GPgh8EF8w3LOuej+9vE6HnxnNRcc1ZfvHts/6HBcHRrScrkKGAGUAU8BhYSf1nfOuYR6e/k2bn5pEV8emsGtZ4zwmWGNWENmixUDN0U+zjkXiMWbCvi/J+cytEcH7v/WEbTwNcMatYbMFhsiabqkNyX9e+8nFjeXdI0kk5Req7yvpCJJ19Qou0PSekn1Ljkj6QZJqyQtl3RKLOJ0zgVr464SvvuXT+jYtiWPXXIk7Vs3pEffBakhf0J/Bx4EHgGqYnVjSVnAScC6KLunAa/XKnsFuB9YWc81hwPnE+7G6w3MkjQksoyNcy4JFRRXcMljsykuq+KZK46mZyefcpwMGpJcKs3sT3G49zTgWuClmoWSJhFeeXlPzXIz+yiyv75rngXMNLMyYI2kVYSnUX8Ys6idcwlTVlnFlBnZrN2+h8e/O84Xo0wiDem0fEXSDyT1ktR17+dgbirpTGCjmeXUKk8DrgNuO8BL9wHW19jeECmLFsMUSdmSsvPy8g7wds65eKmuNn76TA4fr9nBb78xmmMOSd/3Sa7RaEjL5eLIrz+rUWbAwPpOkjQLiLY06U3AjcDJUfbdBkwzs6IDnAUS7aSoz+SY2XRgOkAoFPLndpxrRMyMX766hH8s2Mz1pw3jrMOj/ozoGrGGzBY7oKfxI++A+QJJIwk/4Z8TSSCZwFxJ44DxwGRJU4HOQLWkUjO7v4G33QBk1djOBDYdSPzOueA8+E4uf35/Ld89dgDfm1Dvz7GukaqzW0zStTW+f6PWvl8f6A3NbKGZdTez/mbWn3BCOMLMtpjZ8TXK7wV+vR+JBeBl4HxJrSUNAAYDsw80Vudc4v09ez13/XMZZ47uzc+/dqg/y5Kk6htzOb/G9xtq7Ts1DrHUS9JUSRuAdpI2SLo1Un6mpNsBzGwx8AywhPBCm1f6TDHnksesJVu5/vmFHDcond9+YzQpKZ5YkpXMog83SJpnZmNqf4+2nexCoZBlZ2cHHYZzzdrHudu56LHZDO3Zgb9dfpQ/y5IEJM0xs1C0ffW1XKyO79G2nXPugC3aWMBlj2eT2aUtf/nOOE8sTUB9f4KjJRUSnoHVNvKdyLY/xeSci4ncvCIu+fNsOrRpwYxLx9M1rVXQIbkYqDO5mJm/0s05F1cbd5VwwSMfYwYzLhtP785tgw7JxYi3PZ1zgcjbXcYFj3zM7rJKZk45ikMy2gcdkoshX1bUOZdwu4rLufDRj9lSUMpfvnMkI3r7mySbGk8uzrmEKiyt4KLHZpObt4fpF41lbL+DWk3KNVKeXJxzCbOnrJLv/PkTlmwq5IFvH8HxgzOCDsnFiY+5OOcSoqS8issez2b++l3c/80xTBzeI+iQXBx5y8U5F3elFVVc/kQ2H63Zzu++MZrTRvYKOiQXZ55cnHNxtTexvL86n7snj2bSGF/huDnw5OKci5vSiiq+N2MO763K565zRjF5bGbQIbkE8eTinIuLvYnlnRV53Pn1kZwbytr3Sa7J8AF951zMlZT/ryts6jmjOPdITyzNjScX51xMFZdXculfwoP3d08e7V1hzZQnF+dczOwureC7f/mEOZ/u5J5zR3P2GE8szZUnF+dcTOwqLufix2azeFMhv//mGE4f1TvokFyAPLk45w5aflEZFz46m9XbivjTBWM5yR+QbPY8uTjnDsqmyLL5mwpKeOTiEBOG+JIuzpOLc+4g5OYVceGjsyksqWDGpeM5sr8vQunCPLk45w7I4k0FXPzYbMzgqSlHcVgfXzbf/U+gD1FKukaSSUqvVd5XUpGka2qU3SFpvaSieq7XX1KJpPmRz4PxjN+55uqD1fmc99BHtEpN4ZkrjvbE4r4gsJaLpCzgJGBdlN3TgNdrlb0C3A+s3MelV5vZ4QcdoHMuqtcWbubqmfPp160dT1w6jl6d/NXE7ouCbLlMA64FrGahpElALrC4ZrmZfWRmmxMWnXPuC2Z89ClX/m0uIzM78fcrjvbE4uoUSHKRdCaw0cxyapWnAdcBtx3E5QdImifpHUnH1xPDFEnZkrLz8vIO4nbONX3V1cZd/1zGzS8u4sSh3fnrpePp3K5V0GG5Rixu3WKSZgE9o+y6CbgRODnKvtuAaWZWJOlAbrsZ6Gtm2yWNBV6UNMLMCmsfaGbTgekAoVDIau93zoWVV1Zz7bM5vDh/E98c15dfnjWCFqm+5q2rX9ySi5lNjFYuaSQwAMiJJJBMYK6kccB4YLKkqUBnoFpSqZnd38B7lgFlke9zJK0GhgDZB1kd55qlguIKrvjrHD7M3c41Jw/hyhMGcYA/+LlmJuED+ma2EOi+d1vSWiBkZvnA8TXKbwWKGppYIudkADvMrErSQGAw4fEb59x+Wre9mEv+Mpv1O4r53TdGc44vQOn2Q9K0bSVNlbQBaCdpQyT5IOlMSbdHDpsALJCUAzwLXGFmO4KJ2LnkNefTHUx64H22F5Uz49LxnljcfpOZDzeEQiHLzvaeM+cAnpuzgRueX0jvzm147JIjGZjRPuiQXCMlaY6ZhaLt8yf0nXMAVFUbU99YxkPv5HLMId3447eOoEuazwhzB8aTi3OOwtIKfjxzPm8t28a3x/fl1jNH0NJnhLmD4MnFuWZudV4Rlz+Rzbrtxdx+1gguOrp/0CG5JsCTi3PN2FtLt3L1zPm0apHCXy8bz1EDuwUdkmsiPLk41wxVVRv3zlrBH/69isP6dOShC0P06exLubjY8eTiXDOzc085P5w5j3dX5nNuKJPbzzqMNi1Tgw7LNTGeXJxrRuZ8upOr/jaX/KJy7vz6SM4f1zfokFwT5cnFuWbAzHj43Vym/nM5vTq34dnvH82ozM5Bh+WaME8uzjVxO/aU87O/5/DWsm2cOqInd00eRae2LYMOyzVxnlyca8I+WJXP1U/PZ1dxBb84YziXHNPfF550CeHJxbkmqLyymmmzVvDgO6sZkJ7Gn79zJCN6+6uIXeJ4cnGuiVm1bTdXPz2fRRsLOS+UxS/OHE67Vv5P3SWW/41zromorjae+HAtv3l9GWmtW/DQhWM5ZUS09/U5F3+eXJxrAtbvKOZnz+bwUe4Ovjw0g6mTR9G9Q5ugw3LNmCcX55JYdbXx5Ox1/Oa1paRI3HXOSM4NZfmgvQucJxfnklRuXhHXP7+Q2Wt2cNygdO6aPMqXcHGNhicX55JMRVU1j7y7hntnraBVixSmnjOKb4QyvbXiGhVPLs4lkey1O7jphUUs37qbU0b04JdnHUb3jj624hofTy7OJYEde8q5+41lPDV7Pb07teHhi0KcNLxH0GE5VydPLs41YlXVxt8+/pTfvrmCorJKLj9+AFdPHEJaa/+n6xq3QN9jKukaSSYpvVZ5X0lFkq6JbLeT9KqkZZIWS7qznmveIGmVpOWSTol3HZyLlw9Xb+eMP7zHzS8tZkTvjrz+o+O56WvDPbG4pBDY31JJWcBJwLoou6cBr9cq+62ZvS2pFfCWpNPM7HPHSBoOnA+MAHoDsyQNMbOq2NfAufhYm7+HX7+2lDeXbKVP57bc/60xfG1kLx+wd0klyB+BpgHXAi/VLJQ0CcgF9uwtM7Ni4O3I93JJc4HMKNc8C5hpZmXAGkmrgHHAh/GogHOxlF9Uxv3/XsWTH39Kq9QUfnbKUC49boC/yMslpUCSi6QzgY1mllPzpzFJacB1hFs019RxbmfgDOC+KLv7AB/V2N4QKYt2nSnAFIC+ff2FSS44RWWVPPbeGh56ZzWlldWcG8rixxMH+ywwl9TillwkzQKiLWx0E3AjcHKUfbcB08ysKFoXgKQWwFPA780sN9pto5RZtPjMbDowHSAUCkU9xrl4KimvYsZHa/nTf1azs7iCU0b04GenDGNQ9/ZBh+bcQYtbcjGzidHKJY0EBgB7Wy2ZwFxJ44DxwGRJU4HOQLWkUjO7P3L6dGClmd1bx203AFk1tjOBTQdZFediqqS8iic//pSH/ptL3u4yJgzJ4CcnDeHwrM5Bh+ZczCS8W8zMFgLd925LWguEzCwfOL5G+a1A0d7EIulXQCfgsnou/zLwN0n3EB7QHwzMjnEVnDsghaUVPPnROh55N5fte8o5emA3/vitIxg3oGvQoTkXc0kxp1FSJuHutGWEWzkA95vZI5Hxm5CZ3WJmiyU9AywBKoEr4z1TrKyyitYtfMDV1W1bYSmPvb+WJz/6lN1llUwYksEPTxxEqL8nFdd0ycyHG0KhkGVnZ+/3ebl5RZz+h/f40pAMTh7RgxOH9qBTO383uQtbvKmAx95byys5m6isrua0kb343oSBjMrsHHRozsWEpDlmFoq2LylaLo1Vi5QUvn5EH95cvJXXF20hNUWM69+VicN7MPHQ7vTrlhZ0iC7BKqqq+deSrTz+wVo+XrODdq1SOe/ILC49bgD90/3vg2s+vOXCgbdc9qquNhZsLODNxVuYtXQrK7YWAXBIRhonDuvOCcO6E+rXlVYtAl0QwcXRxl0lPP3JembOXse23WX06dyWS47pz7lHZtGprbdmXdNUX8vFkwsHn1xqW7e9mFlLt/L28m18nLuD8qpq0lqlcuygdL40NIMJgzPI6touZvdzwSirrGLWkm08nb2ed1fmAfClIRlceFQ/vjy0O6kp/kS9a9o8uexDrJNLTXvKKnlvVT7vrMjjneV5bNxVAsCA9DSOH5zOsYPSOWpgN//pNklUVxtz1u3k+bkbeXXBJgpLK+ndqQ2TQ1l8Y2ym/9DgmhVPLvsQz+RSk5mxOq+Id1fm8+7KfD5cvZ2SiipSBIf16cTRh3Tj6IHdCPXvSntfnLDRMDNyNhTw6oJNvLpgM5sKSmnbMpVTD+vJ2WP6cOygdG+luGbJk8s+JCq51FZeWc389bt4f1U+H6zOZ/76XVRUGakp4rA+nRg/oCvj+ncl1L8Lndu1Snh8zVlFVTWfrNnBm0u28ubiLWwqKKVlqpgwOIPTR/fi5OE9fXVi1+x5ctmHoJJLbcXllcz9dBcf5uYze80OctYXUF5VDcCg7u0Z27cLY/t14fC+nRmU0Z4U/2k5prYVlvLOijz+szyP/67MY3dpJa1bpHD84AxOGdGDk4f39KnmztXgyWUfGktyqa20oooFGwrI/nQHc9buZM66newqrgCgQ+sWjMzsxKjMzozO7MTIzE706dzWl2XfD7uKy5m9ZgcfrN7O+6vyWbktPMuve4fWnDA0PMtvwpB02rXyFopz0fhzLkmqTctUxg3o+tnyIGbGmvw9zFu3i7nrdrJgQwGPvpdLRVX4B4RObVsyondHRvTuyLCeHTm0V0cGdW/vU6AJ/96t3V7MvHU7mbtuJ9lrd7Jsy24AWrdIYdyArkwem8lxg9MZ3qujJ2nnDpK3XGi8LZeGKK2oYtmW3SzaWMDiTYUs2ljA8q27Ka8Md6elpogB6WkM6dGeQd07cEhGGodktGdgRlqT/Ym8oqqa3Lw9LNtSyJJNhSzaVMCijYUUlIRbfe1bt+DwrM6MH9CV8QO7MTqrky/h49wB8JZLE9amZSqHZ3X+3Iq6lVXVrMnfw5LNhazYupsVW4tYsqmQ1xdtoebPEj06tqZftzT6dW1HVtd2ZHVtS2aXdvTq1IYeHdvQMrXxtnjMjLyiMtbvKGbdjmLW5O1hdd4eVucVsTqv6LPWXKvUFIb16sBXR/ZiZJ9OHNGvM4O7d/DZXc7FmSeXJqhFagqDe3RgcI8Onysvraji0+3FrM4rIjeviLXbi/l0+x7+syKPvN1lnztWCo89dO/QJvxrx9Z0S2tN17RWdGvfik5tW372ad+mBe1bt6Bty9QD7k4yM8oqqykqq6SotJJdJRXsKi5nV3EF+UVl5BWVkbe7jC0FpWwuKGVzQQmlFdWfnZ8i6Nu1HYdktOdLQzM4tGdHhvXqwMB07xZ0LgieXJqRNi1TGdqzA0N7dvjCvtKKKjbsLGHDzmK2FJSyqaCULQUlbNtdxqaCUnI2FLBjTxnV9fSipghat0ildcsUWrdIoUVKCi1SRYsUfS7pVFUbldXVVFaFE0pZRRVlldVU1nPxVqkppLdvRc9ObRjeuyMTD+1OZpd29O3ajr7d2pHZpa13bTnXiHhycUA48Qzq3r7etyBWVxsFJRVs31NGQUnFZ5+isir2lFWyp6zys2RRWhFOFnuTyF6GkZqSQssUkZqiSCJKpXWLFNJah1tAaa1b0LltS7qktaRzu1akp7WmY9sWPsjuXBLx5OIaLCVFdElrRZc0f6DTOVc/74x2zjkXc55cnHPOxZwnF+ecczHnycU551zMeXJxzjkXc4EmF0nXSDJJ6bXK+0oqknRNZLudpFclLZO0WNKddVyvv6QSSfMjnwcTUQ/nnHOfF9hUZElZwEnAuii7pwGv1yr7rZm9LakV8Jak08ys9jEAq83s8NhG65xzbn8E2XKZBlwLfO6xbEmTgFxg8d4yMys2s7cj38uBuUBmwiJ1zjm3XwJpuUg6E9hoZjk1n7qWlAZcR7hFc00d53YGzgDuq+PyAyTNAwqBn5vZu3VcZwowJbJZJGn5AVQlaOlAftBBJJjXuXlobnVO1vr2q2tH3JKLpFlAzyi7bgJuBE6Osu82YJqZFUVb6kNSC+Ap4Pdmlhvl/M1AXzPbLmks8KKkEWZWWPtAM5sOTG9whRohSdl1LXfdVHmdm4fmVuemWN+4JRczmxitXNJIYACwt9WSCcyVNA4YD0yWNBXoDFRLKjWz+yOnTwdWmtm9ddyzDCiLfJ8jaTUwBEjOl7U451ySSni3mJktBLrv3Za0FgiZWT5wfI3yW4GivYlF0q+ATsBldV1bUgaww8yqJA0EBhMev3HOOZdASfGci6RMwt1pwwm3cuZLuiyy70xJt0cOnQAskJQDPAtcYWY7Agk6MZK6W+8AeZ2bh+ZW5yZXX3/NsXPOuZhLipaLc8655OLJxTnnXMx5ckkikrpK+peklZFfu9RzbKqkeZL+kcgYY60hdZaUJeltSUsjywP9KIhYD4akUyUtl7RK0vVR9kvS7yP7F0g6Iog4Y6kBdf52pK4LJH0gaXQQccbSvupc47gjJVVJmpzI+GLJk0tyuR54y8wGA29FtuvyI2BpQqKKr4bUuRL4qZkdChwFXClpeAJjPCiSUoE/AqcRnrTyzSjxn0Z49uNgwg///imhQcZYA+u8BviSmY0CfkmSD3o3sM57j7sLeCOxEcaWJ5fkchbweOT748CkaAdFZtd9DXgkMWHF1T7rbGabzWxu5Ptuwkm1T6ICjIFxwCozy40sbzSTcL1rOgt4wsI+AjpL6pXoQGNon3U2sw/MbGdk8yOSf8mnhvw5A1wFPAdsS2RwsebJJbn0MLPNEP4PlRrPC9VyL+F126oTFFc8NbTOQHhlbGAM8HH8Q4uZPsD6Gtsb+GJybMgxyWR/63MpX1zMNtnss86S+gBnA0m/ontgqyK76PaxbE5Dzj8d2BZZoeDLMQwtbg62zjWu057wT3xXR1vypxH74lpHtRZ0beAxyaTB9ZF0AuHkclxcI4q/htT5XuC6yIPg8Y8ojjy5NDJ1LZsDIGmrpF5mtjnSJRKt2XwscKakrwJtgI6S/mpmF8Qp5IMWgzojqSXhxPKkmT0fp1DjZQOQVWM7E9h0AMckkwbVR9Iowt27p5nZ9gTFFi8NqXMImBlJLOnAVyVVmtmLCYkwhrxbLLm8DFwc+X4x8FLtA8zsBjPLNLP+wPnAvxtzYmmAfdZZ4X+JjwJLzeyeBMYWK58AgyUNiLyv6HzC9a7pZeCiyKyxo4CCvd2FSWqfdZbUF3geuNDMVgQQY6zts85mNsDM+kf+/T4L/CAZEwt4ckk2dwInSVpJ+LUEdwJI6i3ptUAji5+G1PlY4ELgRP3vLaRfDSbc/WdmlcD/EZ4dtBR4xswWS7pC0hWRw14jvE7eKuBh4AeBBBsjDazzLUA34IHIn2lSL0DbwDo3Gb78i3POuZjzlotzzrmY8+TinHMu5jy5OOecizlPLs4552LOk4tzzrmY8+TinHMu5jy5OOeci7n/B6WbuC/+26WWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(2)\n",
    "plt.plot(d[:,1], d[:,2])\n",
    "plt.ylabel('Energy (eV)')\n",
    "plt.title('Energy')\n",
    "#plt.savefig(\"energy.jpg\")\n",
    "#plt.xlim([0.0,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "NMF (python)",
   "language": "python",
   "name": "nmfenv"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
